scatteredInterpolant in nonlinear system

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 1

[Matt J]

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

[Cristiano Piccinin]

Thanks, Matt for this good idea

Answer 2

[Catalytic]

In the first global declaration, etac is missing from the list.

global Hm pe1 pe2 pe3 etae2 etae3 qmax1 qmax2 qmax3   %<--- add etac

[Cristiano Piccinin]

Thanks, for your comment, you were right and it worked fine.