Algebraic state variables for nlmpc

Hi all,

I am currently developing a nonlinear model predictive control that incorporates algebraic state variables (they are not equipped with derivative equation). Their equations would be in the form of 0 = Fx(k) + Gu(k). I need these state variables to complete the system. I am planning to utilize the Nonlinear Model Predictive Control Toolbox to control the system. However, I am not sure how to add this function in the model function without the derivative. Can anyone kindly explain to me how to incorporate these state variables to the nlmpc code?

Thank you in advance.

Answer 1

[LeoAiE]

The Nonlinear Model Predictive Control Toolbox does support Differential-Algebraic Equation (DAE) systems, and you can incorporate algebraic equations into your model function.

function [dx, y] = myModel(t, x, u, varargin)
    % x: Differential states
    % u: Inputs
    % dx: Derivatives of differential states

    % Define your differential equations
    dx(1,1) = ... % function of x, u and possibly t

    % Define your output equations
    y = ...

    % Incorporate your algebraic equations
    if nargin >= 4
        % varargin{1} contains the algebraic states
        xa = varargin{1};

        % Define algebraic equations
        dx(2,1) = 0 - F*x(1) - G*u; % 0 = Fx + Gu
    end
end

And here is how to setup your nlmpc object with algebraic equations:

nx = 1; % Number of differential states
nu = 1; % Number of inputs
ny = 1; % Number of outputs
nv = 1; % Number of algebraic states
Ts = 0.1; % Sample time
model = @myModel; % Model function

controller = nlmpc(nx, ny, 'Model', model, 'Ts', Ts, 'MV', 1);

% Define number of algebraic states
controller.Model.NumberOfVariables = nx + nu + nv;

% Define Jacobian of the model with respect to algebraic states
controller.Model.Jacobian.OutputVariables = @(t, x, u, varargin) ...
    [0, 0, -G]; % Assuming G is constant w.r.t x, u, and xa

In this way, the nonlinear MPC controller knows that it must solve the algebraic equation for xa at each control interval.

Note that providing the Jacobian can significantly improve the controller's performance. The Jacobian function should return the partial derivatives of the output equations with respect to all state, input, and algebraic variables. Please adjust according to your specific model and function F. Also, ensure the dimensions match your model size.

Remember to use your specific model details in place of this simplified example. Make sure to test this on a smaller scale or with test data to ensure functionality before implementing it into your system. Finally, this is just a starting point. You may need to adjust this code to fit your needs, such as defining prediction and control horizons, constraints, weights, and other specifics of your system and control strategy.

Hope this helps!

[Saskia Putri]

Hi

Thank you very much for providing a clear and detail examples. I have developed the coding following your approach. The ValidateFcn also did not detect any issues.

In summary, I have 2 state vars, 2 control vars, 2 outputs (y=x), 3 measured disturbances, and 2 algebraic states. I have also added the Jacobian Fcn for the state.

However, when I try to run the code, I have the following warning, even though I have set the initial conditions for my state and control variables.

------------------------------------------------------

Start of Error Report

------------------------------------------------------

Error using sqpInterface

Nonlinear constraint function is undefined at initial point. Fmincon cannot continue.

Error in fmincon (line 895)

[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 withHL (line 211)

[uk,nloptions,info] = nlmpcmove(nlobj,xk,lastMV,yref,vk);

------------------------------------------------------

End of Error Report

------------------------------------------------------

Do you have suggestions?

Thank you.

[LeoAiE]

Maybe I can help you if you share your code

[Saskia Putri]

Hi

I have attached the code for your reference. Please kindly let me know if it's unclear.

Thank you.

[LeoAiE]

Here are a few points you should check:

1. Check your initial conditions: You need to ensure that the initial condition you have set for your problem is feasible, i.e., it satisfies all the constraints of your problem. If the initial condition doesn't meet this requirement, the optimizer may not be able to find a feasible direction to move.
2. Check your nonlinear constraints: You should also check your nonlinear constraints. Ensure that they are correctly implemented, and there's no division by zero, square root of a negative number, log of a non-positive number, or anything else that could lead to undefined results.
3. Check your Jacobian function: If you're using a Jacobian, make sure it's correctly defined. The error message suggests that there might be a problem with your function. To verify the correctness of your Jacobian function, you can compare the numerical Jacobian obtained using finite differences to your analytical Jacobian. In MATLAB, you can use the numjac function for numerical Jacobian calculation.

In your code, you are using fmincon indirectly by calling the nlmpcmove function, which in turn calls fmincon. The error seems to occur at the nlmpcmove function, which means you should check your model, constraints, and initial conditions that you have set in your nlmpc object.

In your 'WDSFcn' function, add a check for 'NaN' and 'Inf':

function dx= WDSFcn(x, u, varargin)
    % x: Differential states
    % u: Inputs
    % dx: Derivatives of differential states

    %parameters
    ATank = 153.86; %m^2

    %Pipes
    Dpipe = 0.6096; %m diam of pipe
    apipe = (1/4)*pi*(Dpipe)^2;

    %constant for outflow
    Cd = 1;
    K = Cd*apipe*sqrt(2*9.81)*3600;

    % Define your differential equations
    dx(1,1) = u(1)/ATank -  (K/ATank)*sqrt(x(1)); 
    dx(2,1) = u(2)/ATank -  (K/ATank)*sqrt(x(2));

    % Incorporate your algebraic equations
    if nargin >= 3
        % varargin{1} contains the algebraic states
        p1 = varargin{1};
        p2 = varargin{2};
        

        % Define algebraic equations
        dx(3,1) = 0 - K*sqrt(x(1)) +u(2)-p1+u(3); % 0 = Fx + Gu
        dx(4,1) = 0 - K*sqrt(x(2)) +p2+u(5); % 0 = Fx + Gu
        dx(5,1)= 0 - p2 + p1 + u(4);
    end

    % NEW CODE: Check for Inf and NaN values in dx
    if any(isinf(dx) | isnan(dx))
        error('NaN or Inf found in dx')
    end
end

In your main simulation code, you could add a 'try-catch' block around the 'nlmpcmove' function:

try
    [uk,nloptions,info] = nlmpcmove(nlobj,xk,lastMV,yref,vk);
catch ME
    warning('Error in nlmpcmove: %s', ME.message);
    break;  % Exit the loop
end

This will stop the simulation if 'nlmpcmove' throws an error, and print a warning message.

[Saskia Putri]

Hi

Thank you very much for your advice. I have fixed the code and it is able to run. I had to adjust my dynamic and jacobian fcns. There were some units that are not consistent. Thank you very much, again for the help!

In regards to the algebraic states (p1,p2), may I ask how to extract the data and possibly add them in an equality constraint?

Thank you!

[LeoAiE]

The algebraic states can be accessed within the nlmpcmove function using varargin in the custom ODE function.

However, to use them as part of an equality constraint in your nonlinear MPC problem, they should be included in your nlmpc object's Model.StateFcn and Model.OutputFcn, so that they are part of the optimization problem.

Let's denote your algebraic states as za in the following code.

Here's how you might set up your nlmpc object to include algebraic states and use them as part of an equality constraint:

% Define the numbers of each type of variable in your problem.
nx = 2;  % number of differential states
nu = 5;  % number of inputs
ny = 2;  % number of outputs
nz = 2;  % number of algebraic states

% Create your nlmpc object.
nlobj = nlmpc(nx,ny,nu,nz);

% Specify your prediction model and custom ODE function.
nlobj.Model.StateFcn = @(x,u,za) WDSFcn(x,u,za);
nlobj.Model.OutputFcn = @(x,u) MyOutputFcn(x,u);
nlobj.Model.IsContinuousTime = true;

% Set up equality constraint for your algebraic states.
nlobj.Constraints.p1.Min = nlobj.Constraints.p1.Max = some_value;  % Replace some_value as required
nlobj.Constraints.p2.Min = nlobj.Constraints.p2.Max = some_value;  % Replace some_value as required

% Define the custom ODE function.
function dx = WDSFcn(x, u, za)
    % Your existing code here, but also use za(1) for p1 and za(2) for p2.
end

% Define the custom output function.
function y = MyOutputFcn(x, u)
    % Your existing code here.
end

In the above code, replace some_value with the actual values you need.

This is just a rough guide, and you may need to adjust it based on your specific requirements.

Please if you found the answer helpful consider accepting it!

[Saskia Putri]

I understand. Thank you very much!

[Lars]

Hello,

is it as well possible to define in the state function a system of DAE for the states itself? Something like the following equations:

This system can be solved as specified e.g. as described here: (https://de.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html) But I haven't seen any option to specify a DAE system within the state function for a NLMPC model. Is there any option?

Thanks in advance.