Solver does not recognize problem constraint

[oaldaify]

I hope the discussion title and the code are sufficient for you. The problem is attached for your reference.

I did not add the last two terms of the objective function to the problem yet, but I believe the solver should recognize the constraints on the system variable.

import casadi.*                                
opti = casadi.Opti();       % Create a new optimization problem
opti.solver('ipopt')        % Choose solver

% Parameters
T = 5;                             % continuous-time prediction horizon
N = 20;             %prediction horizon


n = 16;	% State dimension
m = 4;	% Input dimension
x = MX.sym('x',n);
u = MX.sym('u',m);

%% System parameters

%Graph model parameters
N_v = 18;  %No. of vertices
N_e = 21;  %No. of edges
N_in = 2;  %No. of sources
N_out = 2; %No. of sinks

% Define the masses of the components
masses = MX([16; 70; 833; 16; 70; 0; 16; 70; 0; 65; 25; 137; 0; 31; 35; 0]);

%Efficiencies
Eta_batt = 0.90; %Battery 
Eta_c = 0.93;    %Converter 
Eta_g = 0.90;    %Generator 
Eta_r = 0.93;    %ATRU 
Eta_bus = 0.99;  %Bus
Eta_m = 0.90;    %Motor 
Eta_i = 0.93;    %Inverter

As_CP_MD = 3.534;
cp = 0.85;
As_CP = 2.262;
h = 1;

%Heat exchanger
A_s_cool = 60;    %Heat Transfer Area Coolant
A_s_air = 103;    %Heat Transfer Area Air

rho_air = 1;
V = 1;
A_c = 0.1;
m_fan = 5;

%Battery
C_batt = 250;        %Capacity
C_rate = 1/5;        %Discharging rate

% Define the capacitances of the components
capacities = [cp*masses(1), cp*masses(2), C_batt, ...
              cp*masses(4), cp*masses(5), 1, ...
              cp*masses(7), cp*masses(8), 1, ...
              cp*masses(10), cp*masses(11), cp*masses(12), ...
              1, cp*masses(15), cp*masses(14), 1];

%Weighting matrix
Wm = diag(1e7* [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0]);
Wu = [1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0]; %For mass flow rate  in the TMS
Wu_dot = 1e9 * Wu;

% Create the capacitance matrix
P_prop = 1500;
C = diag(capacities);

M = [-1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
     0 -1 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0;
     1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
     0 0 0 -1 -1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
     0 0 -1 0 0 1 1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 -1 -1 0 0 0 0 1 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 -1 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -1 1 0 0 0 0 0;
     0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 0 0 1;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 1 0 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1];

M_bar = M(1:N_v-N_out, :);

D = [0 0; 0 0; 0 0; 0 0; 0 0; 0 0; 0 0; 0 0;
     0 1; 0 0; 0 0; 0 0; 0 0; 0 0; 0 0; 1 0;];
 
P = [u(1)*(1-Eta_batt*Eta_c)/(Eta_batt*Eta_c);
    h*As_CP*(x(1)-x(2));
    u(1);
    u(3)*cp*x(2);
    h*As_CP*(x(4)-x(5));
    (Eta_g*Eta_r*u(2)+Eta_c*u(1))*(1-Eta_bus)/(Eta_m*Eta_i);
    (Eta_g*Eta_r*u(2)+Eta_c*u(1))/(Eta_m*Eta_i);
    u(3)*cp*x(10);
    u(3)*cp*x(5);
    h*As_CP*(x(7)-x(8));
    u(2)*(1-Eta_g*Eta_r);
    u(2)*Eta_g*Eta_r;
    (Eta_g*Eta_r*u(2)+Eta_c*u(1))*(1-Eta_m*Eta_i)/(Eta_m*Eta_i);
    h*As_CP_MD*(x(11)-x(12));
    u(3)*cp*x(7);
    u(3)*cp*(12);
    u(3)*cp*x(14);
    h*A_s_cool*(x(14)-x(15));
    h*A_s_air*(x(15)-x(16));
    u(4)*cp*x(15);
    (Eta_g*Eta_r*u(2)+Eta_c*u(1))];

P_in = [u(4)*cp*x(16); u(2)];
P_track = 500;
x_dot = inv(C) * (D*P_in - M_bar*P);

f = Function('f',{x,u},{x_dot},{'x','u'},{'x_dot'});
%Construct integration to discretize the system
G = struct;
G.tf = T/N; %integratin period
G.simplify = true;

%Construct DAE
dae = struct;
dae.x = x; %define the states
dae.p = u; %define the control parameter
dae.ode = f(x,u); %define the RH side of ode
intg = integrator('intg','rk',dae, G); %Integration command

masses = opti.variable(n,N);
u = opti.variable(m, N);
x = opti.variable(n, N+1);
u_dot = opti.variable(m,N);
% Constraints
% dynamic constraint
for k=1:N    
   opti.subject_to(x(:,k+1) == f(x(:,k),u(:,k)));
end

%Fluid temperature constraint
opti.subject_to(20< x(2,:) <= 80)
opti.subject_to(20<x(5,:) <= 80)
opti.subject_to(20<x(8,:) <= 80)
opti.subject_to(20<x(10,:) <= 80)
opti.subject_to(20<x(12,:) <= 80)
opti.subject_to(20<x(14,:) <= 80)

%CP wall temperature constraint
opti.subject_to(20<x(1,:) <= 50)
opti.subject_to(20<x(4,:) <= 90)
opti.subject_to(20<x(7,:) <= 90)
opti.subject_to(20<x(11,:) <= 90)

opti.subject_to(0 <u(2,:)<1.4*P_track)
opti.subject_to(0.1 <= x(3,:) <= 1)

opti.subject_to(rho_air * V* A_c <= u_dot(4,:) <= (rho_air * V* A_c) +m_fan)
opti.subject_to(0 <= u_dot(3,:) <= 8)

opti.subject_to(u(1,:) - C_batt*C_rate == 0)

opti.subject_to((Eta_g*Eta_r*u(2,:)+Eta_c*u(1,:) - (P_prop/(Eta_bus*Eta_m*Eta_i))) == 0)

% cost function
J = 0;
Ji = 0;
Ju = 0;


for i=1:N
    J = masses(:,i)'*Wm*masses(:,i);
     Ji = Ji + u(:,i)'*Wu*u(:,i);

     J_m = J+Ji;
end

opti.minimize(J_m);

% Initial guess
u_init = ones(m,N);            % Initial guess for input sequence: all 0
opti.set_initial(u,u_init);
x_init = ones(n,N+1);
opti.set_initial(x,x_init);             % Initial guess for state sequence via simulation

% solve the OCP
sol = opti.solve();
x_OL = sol.value(x);
u_OL = sol.value(u);
m_OL = sol.value(masses)

[svigerske]

Maybe that is something for the CASADI people to look at first.

[JaffaSweet]

good thanks