MATLAB Answers

scatteredInterpolant in nonlinear system

조회 수: 2(최근 30일)
Cristiano Piccinin
Cristiano Piccinin 2019년 6월 8일
편집: Matt J 2019년 6월 10일
Hi,
I'm trying to implement solution of a nonlinear system, in which i'd like to use a scatteredInterpolant to calculate some values. 9 equations.
I have 3 vectors containing the data in which to construct the interpolant: Q, H and P_idr_gr1 , that is power of turbine vs. discharge and head.
The interpolant used from the command window works just fine.
When i try to implement in the nonlinear system i get the errors:
Index in position 1 is invalid. Array indices must be positive integers or logical values. (or: Index in position 1 exceeds array bounds, depending on x0)
Error in nlsystem>mysystem (line 55)
S(7) = x(8)*x(4)-pe1*1000/gamma/etac(x(8),x(4));
Error in fsolve (line 242)
fuser = feval(funfcn{3},x,varargin{:});
Error in nlsystem (line 41)
[x1,fval1,exitflag1] = fsolve(@mysystem,x0,options);
Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.
I'd really appreciate to understand what i'm, doing wrong !
Thanks,
Cristiano
% my nonlinear system using scatteredInterpolant
clear all
clc;
global Hm pe1 pe2 pe3 etae2 etae3 qmax1 qmax2 qmax3
%% dati iniziali
load '190607 collinare USBR fig 35 rendimento turbina.mat'
Hm = 485.90;
pe1 = 10.0;
pe2 = 1.40;
pe3 = 1.40;
etae_2 = [-0.0131 0.042 0.0138 0.7827]; etae2 = polyval(etae_2, pe2);
etae_3 = [-0.0131 0.042 0.0138 0.7827]; etae3 = polyval(etae_3, pe3);
% here i create surface of hill chart
Q100 = 27.9;
H100 = 81;
Q = Q100.*cUSBRfig35.rel_q;
H = H100.*cUSBRfig35.rel_h;
P_idr_gr1 = 9.8045*Q.*H;
% remove generator losses and calculate efficiency eta
P_erog_gr1 = P_idr_gr1.*cUSBRfig35.eta_t - (0.0089.*P_idr_gr1+291.77);
eta_erog_gr1 = P_erog_gr1./P_idr_gr1;
% calculate interpolated surface
etac = scatteredInterpolant(Q,H,eta_erog_gr1,'linear');
% some initial condition for fsolve
x0 = [ Hm Hm Hm (Hm-387.29) (Hm-387.29) 387.29 0.0 0.0 0.0 ];
%% lancio solutore
options = optimoptions('fsolve','Display','iter', 'FunctionTolerance', 1e-6);
[x1,fval1,exitflag1] = fsolve(@mysystem,x0,options);
function S = mysystem(x)
global Hm pe1 pe2 pe3 etae2 etae3 etac
gamma = 9.8045;
S(1) = Hm -x(1) -0.018588*x(7)^2;
S(2) = x(1)-x(2)-0.000541*x(8)^2;
S(3) = x(1)-x(3)-0.007589*x(9)^2;
S(4) = x(4)-x(2)+x(6);
S(5) = x(5)-x(3)+x(6);
S(6) = x(6) - (0.0003*x(7)^3 - 0.0183*x(7)^2 + 0.4374*x(7) + 387.29);
S(7) = x(8)*x(4)-pe1*1000/gamma/etac(x(8),x(4));
S(8) = x(9)*x(5)-(pe2*1000/gamma/etae2 + pe3*1000/gamma/etae3);
S(9) = x(7)-x(8)-x(9);
end

채택된 답변

Matt J
Matt J 2019년 6월 8일
편집: Matt J 2019년 6월 10일
Try this version, which uses nested functions instead.
function nlsystem
% my nonlinear system using scatteredInterpolant
%% dati iniziali
Lstruct=load('190607 collinare USBR fig 35 rendimento turbina.mat');
cUSBRfig35 = Lstruct.cUSBRfig35;
Hm = 485.90;
pe1 = 10.0;
pe2 = 1.40;
pe3 = 1.40;
etae_2 = [-0.0131 0.042 0.0138 0.7827]; etae2 = polyval(etae_2, pe2);
etae_3 = [-0.0131 0.042 0.0138 0.7827]; etae3 = polyval(etae_3, pe3);
% here i create surface of hill chart
Q100 = 27.9;
H100 = 81;
Q = Q100.*cUSBRfig35.rel_q;
H = H100.*cUSBRfig35.rel_h;
P_idr_gr1 = 9.8045*Q.*H;
% remove generator losses and calculate efficiency eta
P_erog_gr1 = P_idr_gr1.*cUSBRfig35.eta_t - (0.0089.*P_idr_gr1+291.77);
eta_erog_gr1 = P_erog_gr1./P_idr_gr1;
% calculate interpolated surface
etac = scatteredInterpolant(Q,H,eta_erog_gr1,'linear');
% some initial condition for fsolve
x0 = [ Hm Hm Hm (Hm-387.29) (Hm-387.29) 387.29 0.0 0.0 0.0 ];
%% lancio solutore
options = optimoptions('fsolve','Display','iter', 'FunctionTolerance', 1e-6);
[x1,fval1,exitflag1] = fsolve(@mysystem,x0,options);
function S = mysystem(x)
gamma = 9.8045;
S(1) = Hm -x(1) -0.018588*x(7)^2;
S(2) = x(1)-x(2)-0.000541*x(8)^2;
S(3) = x(1)-x(3)-0.007589*x(9)^2;
S(4) = x(4)-x(2)+x(6);
S(5) = x(5)-x(3)+x(6);
S(6) = x(6) - (0.0003*x(7)^3 - 0.0183*x(7)^2 + 0.4374*x(7) + 387.29);
S(7) = x(8)*x(4)-pe1*1000/gamma/etac(x(8),x(4));
S(8) = x(9)*x(5)-(pe2*1000/gamma/etae2 + pe3*1000/gamma/etae3);
S(9) = x(7)-x(8)-x(9);
end
end
  댓글 수: 1
Cristiano Piccinin
Cristiano Piccinin 2019년 6월 8일
Thanks, Matt for this good idea

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Catalytic
Catalytic 2019년 6월 8일
In the first global declaration, etac is missing from the list.
global Hm pe1 pe2 pe3 etae2 etae3 qmax1 qmax2 qmax3 %<--- add etac
  댓글 수: 1
Cristiano Piccinin
Cristiano Piccinin 2019년 6월 10일
Thanks, for your comment, you were right and it worked fine.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by