How to solve the error "Error using sqpInterface Nonlinear constraint function is undefined at initial point. Fmincon cannot continue." Error occurred when calling NLP solver

Please I need help to resolve the "Error using sqpInterface Nonlinear constraint function is undefined at initial point. Fmincon cannot continue." Error occurred when calling NLP solver"

I follow the example "Swing-up control of a pendulum using nonlinear model predictive control" to design nonlinear model predictive for a system. I followed the example and adapted it to my system. I created all the files as given in the example. However, I am getting the error above each time i run the code. Nontheless, the exampe runs fine. I am wondering what could be causing the problem.

Here is the code:

nx=13;
ny=1;
nu=1;
nlobj=nlmpc(nx,ny,nu);
Ts=0.1;
nlobj.Ts=Ts;
nlobj.PredictionHorizon=10;
nlobj.ControlHorizon=5;
nlobj.Model.StateFcn="Holos_microreactorDT0";
nlobj.Model.IsContinuousTime = false;
nlobj.Model.NumberOfParameters = 1;
nlobj.Model.OutputFcn = 'Holos_microreactorOutputFcn';
nlobj.Jacobian.OutputFcn = @(x,u,Ts) [1 0 0 0 0 0 0 0 0 0 0 0 0];
nlobj.Weights.OutputVariables = 3;
nlobj.Weights.ManipulatedVariablesRate = 0.1;
nlobj.OV.Min = -10;
nlobj.OV.Max = 150;
nlobj.MV.Min = -0.145;
nlobj.MV.Max = 0.145;
x0 = [0.6;0.1;0.1;0.1;0.1;0.1;0;0;0;0;0;0;0];
u0 = 0.1;
validateFcns(nlobj,x0,u0,[],{Ts});
EKF = extendedKalmanFilter(@Holos_microreactorStateFcn, @Holos_microreactorMeasurementFcn);
x = [1;0;0;0;0;0;0;0;0;0;0;0;0];
y = x(1);
EKF.State = x;
mv = 0;
yref = 0.8;
nloptions = nlmpcmoveopt;
nloptions.Parameters = {Ts};
Duration = 20;
hbar = waitbar(0,'Simulation Progress');
xHistory = x;
for ct = 1:(20/Ts)
    % Set references
    %if ct*Ts<10
    yref;
    %else
    %yref = yref2;
    %end
    % Correct previous prediction using current measurement.
    xk = correct(EKF, y);
    % Compute optimal control moves.
    %z=[0 0 0 0 0 0 0 0 0 0 0 0 0 ]'
    [mv,nloptions,info] = nlmpcmove(nlobj,xk,mv,yref,[],nloptions);
    
    % Predict prediction model states for the next iteration.
    predict(EKF, [mv; Ts]);
    
    % Implement first optimal control move and update plant states.
    x = Holos_microreactorDT0(x,mv,Ts);
    
    % Generate sensor data with some white noise.
    y = x(1) + randn(1)*0.01; 
    
    % Save plant states for display.
    xHistory = [xHistory x]; %#ok<*AGROW>
    waitbar(ct*Ts/20,hbar);
end
close(hbar)

My system has 13 states, single output and single input as given above.

I would appreciate your help, thank you

Answer 1

[Emmanouil Tzorakoleftherakis]

The dynamics/state function are turned into constraints internally when creating a NLP for fmincon. You don't provide all the functions required to run the code so I am not 100% sure but I would guess the error comes from the state function itself. Maybe put a break point and check the output of this function for the specific (x,u) pair that creates the problem.

[Kamal]

@Emmanouil Tzorakoleftherakis Thank you for the insight. I have attached the files below. Yeah I believe the challenge is implict but here are the functions I generated based on my system in the attached files below . I used the same format in the matlab example of nonlinear model predictive of inverted pendulum. I would be glad if this error is spoted and my code runs. Please find attached the functions below. thank you

nx=13;
ny=1;
nu=1;
nlobj=nlmpc(nx,ny,nu);
Ts=0.1;
nlobj.Ts=Ts;
nlobj.PredictionHorizon=10;
nlobj.ControlHorizon=5;
nlobj.Model.StateFcn="Holos_microreactorDT0";
nlobj.Model.IsContinuousTime = false;
nlobj.Model.NumberOfParameters = 1;
nlobj.Model.OutputFcn = 'Holos_microreactorOutputFcn';
nlobj.Jacobian.OutputFcn = @(x,u,Ts) [1 0 0 0 0 0 0 0 0 0 0 0 0];
nlobj.Weights.OutputVariables = 3;
nlobj.Weights.ManipulatedVariablesRate = 0.1;
nlobj.OV.Min = -10;
nlobj.OV.Max = 150;
nlobj.MV.Min = -0.145;
nlobj.MV.Max = 0.145;


x0 = [0.6;0.1;0.1;0.1;0.1;0.1;0;0;0;0;0;0;0];
u0 = 0.1;
validateFcns(nlobj,x0,u0,[],{Ts});

Model.StateFcn is OK.
Model.OutputFcn is OK.
Jacobian.OutputFcn is OK.
Analysis of user-provided model, cost, and constraint functions complete.

EKF = extendedKalmanFilter(@Holos_microreactorStateFcn, @Holos_microreactorMeasurementFcn);
x = [1;0;0;0;0;0;0;0;0;0;0;0;0];
y = x(1);
EKF.State = x;
mv = 0;

yref = 0.8;

nloptions = nlmpcmoveopt;
nloptions.Parameters = {Ts};

Duration = 20;
%hbar = waitbar(0,'Simulation Progress');
xHistory = x;
for ct = 1:(20/Ts)

    % Set references
    %if ct*Ts<10
        yref;
    %else
        %yref = yref2;
    %end

    % Correct previous prediction using current measurement.
    xk = correct(EKF, y);

    % Compute optimal control moves.
    %z=[0 0 0 0 0 0 0 0 0 0 0 0 0 ]'
    [mv,nloptions,info] = nlmpcmove(nlobj,xk,mv,yref,[],nloptions);
    
    % Predict prediction model states for the next iteration.
    predict(EKF, [mv; Ts]);
    
    % Implement first optimal control move and update plant states.
    x = Holos_microreactorDT0(x,mv,Ts);
    
    % Generate sensor data with some white noise.
    y = x(1) + randn(1)*0.01; 
    
    % Save plant states for display.
    xHistory = [xHistory x]; %#ok<*AGROW>
    %waitbar(ct*Ts/20,hbar);
end

------------------------------------------------------
                Start of Error Report                
------------------------------------------------------
Error using sqpInterface
Nonlinear constraint function is undefined at initial point. Fmincon cannot continue.

Error in fmincon (line 900)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = sqpInterface(funfcn,X,full(A),full(B),full(Aeq),full(Beq), ...

Error in nlmpc/nlmpcmove (line 174)
                [z, cost, ExitFlag, Out] = fmincon(CostFcn, z0, A, B, [], [], zLB, zUB, ConFcn, fminconOpt);

Error in LiveEditorEvaluationHelperEeditorId (line 54)
    [mv,nloptions,info] = nlmpcmove(nlobj,xk,mv,yref,[],nloptions);

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON
------------------------------------------------------
                End of Error Report                
------------------------------------------------------

Error using nlmpc/nlmpcmove
Error occurred when calling NLP solver "fmincon". See the error report displayed above.

close(hbar)


figure

subplot(2,2,1)
plot(0:Ts:Duration,xHistory(1,:))
xlabel('time')
ylabel('z')
title('cart position')

subplot(2,2,2)
plot(0:Ts:Duration,xHistory(2,:))
xlabel('time')
ylabel('zdot')
title('cart velocity')

subplot(2,2,3)
plot(0:Ts:Duration,xHistory(3,:))
xlabel('time')
ylabel('theta')
title('pendulum angle')

subplot(2,2,4)
plot(0:Ts:Duration,xHistory(4,:))
xlabel('time')
ylabel('thetadot')
title('pendulum velocity')

function dx = Holos_microreactorCT0(x,u)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% States
% States
n_r = x(1);
Cr1 = x(2); Cr2 = x(3); Cr3 = x(4); Cr4 = x(5);Cr5=x(6);Cr6=x(7);
X = x(8); I = x(9); Tf = x(10); Tm = x(11); Tc = x(12); Rho_r = x(13);


%  Control Inputs
% Zr1 = 0;
% Zr2 = 0;% Output
dx = zeros(13,1) ;

%% Parameters
%D = 0.16;
%v = 2200;
Sig_x=2.65e-22;
yi=0.061;
yx=0.002;
%gama_d=-1.2;
lamda_x = 2.09e-5;
lamda_I = 2.87e-5;
%Lemda_H=0.5;
Sum_f = 0.3358;
%nv=2.2e+3; % average velocisty of neutrons (thermal)
%nu=2.43; % the total number of neutrons liberated per rx
%Sigf=1/(gen*nv*nu); % what is this 




l=1.68e-3;
beta = 0.0065;
beta_1 = 2.267510e-4;
beta_2 = 1.22023e-3;
beta_3 = 1.193061e-3;
beta_4=2.785813e-3;
beta_5=1.257395e-3;
beta_6=5.226188e-4;
Lamda_1 = 0.0124;
Lamda_2 = 0.0369;
Lamda_3 = 0.632;
Lamda_4=3.044508e-1;
Lamda_5=8.539933e-1;
Lamda_6=2.868585;

%Gr8   =-2350E-5/180*2; % All drums
Gr4   =-1450E-5/180*2; % half
%Gr2   =-660E-5/180*2; % 1/4
Gr1   =-250E-5/180*2; % one
%Gr = 0.01;

cp_f=977;
cp_m=1697;
cp_c=5188.6;

M_f=2002;
M_m=11573;
M_c=500;

mu_f=M_f*cp_f;
mu_m=M_m*cp_m;
mu_c=M_c*cp_c;

f_f = 0.96;
P_0 = 20;

Tf0=1105;
Tm0=1087;
T_in=864;
T_out=1106;
%T_c_in=864;
Tc0=(T_in+T_out)/2;
K_fm=f_f*P_0/(Tf0-Tm0);
K_mc=P_0/(Tm0-Tc0);


M_dot=1.75E+01;

alpha_f=-2.875e-5;
alpha_m=-3.696e-5;
alpha_c=0.0;

%% Dynamics
%%kinetics equations with three delayed neutron groups
%rho=Rho_r+gama_d*Q;
rho = Rho_r+alpha_f*(Tf-30)+alpha_c*(Tc-30)+alpha_m*(Tm-30)-Sig_x*(X-1)/Sum_f;
dx(1) = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
dx(2) = Lamda_1*n_r-Lamda_1*Cr1;
dx(3) = Lamda_2*n_r-Lamda_2*Cr2;
dx(4) = Lamda_3*n_r-Lamda_3*Cr3;
dx(5) = Lamda_4*n_r-Lamda_4*Cr3;
dx(6) = Lamda_5*n_r-Lamda_5*Cr3;
dx(7) = Lamda_6*n_r-Lamda_6*Cr3;
%% thermal–hydraulics model of the reactor core
%dQ=-Lemda_H*Q+n_r-n_r0; % n_r0 is the initial power
dx(10)=f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
dx(11)=(1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
dx(12)=K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in);
%% changes of the xenon concentration
G = 3.2e-11;
V = 400*200;% i need to confirm the volume to be used 
Pi = P_0/(G*Sum_f*V);
dx(8) = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
dx(9) = yi*Sum_f*Pi-lamda_I*I;
%dRho_r = Gr*Zr;
dx(13) = (Gr1+Gr4)*u;
end
function xk1 = Holos_microreactorDT0(xk, uk, Ts)
% Repeat application of Euler method sampled at Ts/M.
M = 10;
delta = Ts/M;
xk1 = xk;
for ct=1:M
    xk1 = xk1 + delta*Holos_microreactorCT0(xk1,uk);
end
end
function y = Holos_microreactorMeasurementFcn(xk)
y = xk(1);
end
function y= Holos_microreactorOutputFcn(x, u,params)
y = x(1);
end
function xk1 = Holos_microreactorStateFcn(xk, u)
uk = u(1);
Ts = u(2);
xk1 = Holos_microreactorDT0(xk, uk, Ts);
end

[Emmanouil Tzorakoleftherakis]

Thanks for the files. I did take a look and I am positive the error is caused by the Holos_microreactorCT0 function. Essentially, after the first iteration in your for loop, the state values are the following:

x =

   1.0e+52 *

   -7.9467
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
   -0.0000
    0.0000
    0.0000
   -0.0000
   -0.0000
   -0.0000
         0
         

So when you place such big values in a cost function that is a quadratic, the cost values will become inf and the optimization will complain. You should consider scaling your dynamics to avoid this problem.

Hope this helps

[Kamal]

@Emmanouil Tzorakoleftherakis Thank you for the time and th effort. I will definitely in the direction. You have given me a headway

[Emmanouil Tzorakoleftherakis]

Np. If this answers your question, don't forget to accept the answer. It helps the community if they encounter similar situations.

[Torsten]

@Emmanouil Tzorakoleftherakis

Do you know how to call the internally generated Objective and Constraint functions for "fmincon" with the initial value vector x0 directly ?

[Kamal]

@Emmanouil Tzorakoleftherakis I could not agree more with you. I am trying to rescale my dynamics. I have tried to use rescale function in matlab but error the rescale function seems not working. Could you please point out how to scale the dyanmics in the holomicroreactorCT0? thank you

[Emmanouil Tzorakoleftherakis]

@Torsten this functionality is not exposed. This is part of the value of Model Predictive Control Toolbox - taking the MPC formulation and turning it into an optimization problem. That said, if you evaluate the state function a few times (state function becomes a constraint under the hood), you expect to see some non numerical value here. Also, x0 is not the problem actually here. The problem is x1 (the next time step). Remember MPC solves an optimization problem in receding horizon. With that in mind, the using x and u after the first iteration of the loop before, calling the state function you get

>> Holos_microreactorDT0(x, mv, Ts)

ans =

   1.0e+04 *

       NaN
       NaN
       NaN
       NaN
       NaN
       NaN
       NaN
    0.3125
    9.5312
       NaN
       NaN
       NaN
         0

which explains the error.

@Kamal the problem might be more challenging that simply using rescale. Each state has its own value span, so it would be beneficial to introduce scaling factors for each state and normalize their values in a reasonable set. See the main idea here and here (note this is for linear MPC but you are using the same quadratic cost in your problem). For examples on how to use scale factor in nonlinear MPC see here. Essentially, each state will have its own scale factor which will show up in the const, in the constraints.

Also, please make sure that your dynamics are correct. I cannot comment on that, as this is not my area of expertise.

Hope this helps

[Kamal]

@Emmanouil Tzorakoleftherakis This is great. Thank you for all these. I will explore the links and i belive I should be able to solve the issue now

[Kamal]

@Emmanouil Tzorakoleftherakis Thank you so much. My model is now working excellently. Your advice really helped. I got further help from another staff who visited my university. We got it working but extremely slow. For instance, 10 seconds simulation used to take about 1 hour. I later discovered an error that was causing the the output to be exploding. Immediately I fixed error, the simulation is very fast now and giving accurate and expected results.

However, I dont know why the simulink version is not working? I created bus paramter and other necessary requirements but not working.

2. Could you please advice on how to test for disturbance rejection using the Matlab script of the nonlinear MPC model that is working now? I am trying to pass a sinusoidal disturbance into the system but I dont seem to know how?

[Kamal]

here is the correct model and some of the outputs generated. I have attached the correct model as well as a good simulation result i got. I have not been able to implement the simulink version. Thnak you

[Emmanouil Tzorakoleftherakis]

Np. Can you please post a separate question for this? Otherwise it will get lost in this thread.

Thanks!

Answer 2

[Torsten]

Before calling "fmincon", call the constraint function with your initial vector. My guess is that it returns NaN or Inf for some equality/inequality constraints.

[Kamal]

@Torsten Thank you fro the insight. I didn't write function for constraints since I have passed the constraints on the MV and MO in the body of the code. see below:

nlobj.OV.Min = -10;

nlobj.OV.Max = 150;

nlobj.MV.Min = -0.145;

nlobj.MV.Max = 0.145;

Do I still need to write a function for the constraints as a separate file? Although no constraint function is written in the example I followed? There are like five other files that the above script calls each time it is run.

Thank you and your advice will be appreciated.

[Torsten]

It seems your problem is internally transformed and handed to "fmincon". I don't know how to find the reason for failure if the objective and constraint functions that were internally generated cannot be made accessible from your code. Maybe it's possible, but I have no experience with the MPC Toolbox.

[Kamal]

@Torsten Thank you. Yeah the problem is implict because there is no fmicon in my code. However, I believe the nmpc function needs to call its own execution.