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

Question

Kamal 2024-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073406-how-to-solve-the-error-error-using-sqpinterface-nonlinear-constraint-function-is-undefined-at-initi

评论： Emmanouil Tzorakoleftherakis 2024-2-28

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Emmanouil Tzorakoleftherakis 2024-1-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073406-how-to-solve-the-error-error-using-sqpinterface-nonlinear-constraint-function-is-undefined-at-initi#answer_1396741

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.

11 个评论
显示 9更早的评论隐藏 9更早的评论

Kamal 2024-1-24

编辑：Torsten 2024-1-24

在 MATLAB Online 中打开

@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 2024-1-26

编辑：Emmanouil Tzorakoleftherakis 2024-1-26

在 MATLAB Online 中打开

@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 2024-2-27

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 2024-2-28

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

Thanks!

请先登录，再进行评论。

Answer 2

Torsten 2024-1-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2073406-how-to-solve-the-error-error-using-sqpinterface-nonlinear-constraint-function-is-undefined-at-initi#answer_1395716

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2024-1-23

编辑：Torsten 2024-1-23

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 2024-1-23

@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.

请先登录，再进行评论。