scatteredInterpolant in nonlinear system

3 views (last 30 days)
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

Accepted Answer

Matt J
Matt J on 8 Jun 2019
Edited: Matt J on 10 Jun 2019
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

More Answers (1)

Catalytic
Catalytic on 8 Jun 2019
In the first global declaration, etac is missing from the list.
global Hm pe1 pe2 pe3 etae2 etae3 qmax1 qmax2 qmax3 %<--- add etac
  1 Comment
Cristiano Piccinin
Cristiano Piccinin on 10 Jun 2019
Thanks, for your comment, you were right and it worked fine.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!