scatteredInterpolant in nonlinear system

Asked by Cristiano Piccinin

Cristiano Piccinin (view profile)

on 8 Jun 2019
Latest activity Edited by Matt J

Matt J (view profile)

on 10 Jun 2019
Accepted Answer by Matt J

Matt J (view profile)

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

Answer by Matt J

on 8 Jun 2019
Edited by Matt J

Matt J (view profile)

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

1 Comment

Cristiano Piccinin

Cristiano Piccinin (view profile)

on 8 Jun 2019
Thanks, Matt for this good idea

Catalytic (view profile)

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 (view profile)

on 10 Jun 2019
Thanks, for your comment, you were right and it worked fine.