Cannot create solver

[OJ adukwu]

Good day house,

I have applied this code to a simple 2 states ODE CSTR system and it worked well but I am trying to extend it to a 3 states DAE gas lift  system but I kept getting the error message below;

Error using casadi.nlpsol (line 944)

.../casadi/core/nlpsol.cpp:120: Cannot create 'solver' since [wiv, wrg, wpg, wro, wpo] are free.

These "wiv, wrg, wpg, wro, wpo" defines my states but the solver can not see them. I tried to create another sets of contraints g3 to g7 that defines "wiv, wrg, wpg, wro, wpo"  so that the solver will will see them but I still get the same message.

What is the problem and what is the way out?

I am still very new to programing and CASADI

Thanks in advance.

Below is the code.

clear all

close all

clc

addpath('C:\Users\HP\Documents\casadi')

import casadi.*

% CasADi v3.1.1

T = 60; %[s]x 

 

N = 30;                                % Control horizon 

Q = 1e1*diag([1 1 1]);                 % State weights

Qu = diag([1]*1e1);                  % Input target weights

Ru  = diag([1e-1]);  

Lbh = 500; Lw=1500;  La = 1500;  Hw=1500;  Dw=0.121;  Da=0.189;

 Aw = 0.25*pi*Dw^2;  Aa=0.25*pi*(Da^2-Dw^2);  Va=Aa*La;  Abh = Aw;

rho_o = 800; Ta = 301; Tw = 305;  Mw = 0.028; R = 8.314;

 Civ = 1*1e-4; Cpc = 2e-3; Cr = 2.623*10^-4;

 g = 9.8; wgl = .4; GOR = 0.001;

 Pr = 150*1e5; Ps= 20*1e5; 

 

 

%  wgl = .5;

 

u_min  = [0]';

u_max  = [1]';

du_max = [.05]';

x1max=2000; x1min=1000;

x2max=400; x2min=0;

x3max=10000; x3min=4000;

 

 

%% Initial conditions and set-points

 

x10=1600; x20=200;x30=7000;

delta0=[0;0;0];

u_max = 1; u_min =0;

Valve=SX.sym('Valve'); Pa= SX.sym('Pa'); rho_a = SX.sym('rho_a'); Pwh = SX.sym('Pwh');

rho_w = SX.sym('rho_w'); wpc = SX.sym('wpc'); Pw = SX.sym('Pw'); Pbh = SX.sym('Pbh');

wiv = SX.sym('wiv'); wpg = SX.sym('wpg'); wpo = SX.sym('wpo'); wro = SX.sym('wro');

wrg = SX.sym('wrg');

 

x1 = SX.sym('x1'); x2 = SX.sym('x2');x3 = SX.sym('x3');

states = [x1;x2;x3]; n_states = length(states);

 

u = SX.sym('u');

controls = [u]; n_controls = length(controls);

rhs = [ (wgl-wiv);(wiv+wrg-wpg);(wro-wpo)]; % system r.h.s

rhs1 = [ Civ*sqrt(rho_a*(max (0.0,(Pa-Pw)))); Cr*sqrt(rho_o*(max(0.0,(Pr-Pbh))));(GOR*wro);((x2./max(0.0,(x2+x3))).*wpc);((x3./max(0.0,(x2+x3))).*wpc);];

 

f = Function('f',{states,controls},{rhs}); 

% f1 = Function('f1',{input1},{rhs1});

U = SX.sym('U',n_controls,N); % Decision variables (controls)

P = SX.sym('P',2*n_states + 2*n_controls); % parameters (which include the initial state of the gas lift, the reference state, initial input and desired input)

 

X = SX.sym('X',n_states,(N+1));% A vector that represents the states over the optimization problem.

 

obj = 0; % Objective function

g1 = [];  % constraints vector

 

st  = X(:,1); % initial state

g1 = [g1;st-P(1:3)]; % initial condition constraints;              But could they have subtracted?

 g2 = [];

% g3=[];

% g4=[];

% g5=[];

% g6=[];

% g7=[];

u_km1 = P(2*n_states+1:2*n_states+n_controls);

udes = P(2*n_states+2*n_controls:end);

 udes = 0.95;

x1=x10;

x2=x20;

x3=x30;

for k = 1:N

    

    st = X(:,k);  con = U(:,k);%dcon=con-U(:,k-1); 

    

    cons=udes-U(:,N);

    delta_u = con - u_km1;

 %obj = obj+(st-P(3:4))'*Q*(st-P(3:4)) %+cons*Qu*cons; % calculate obj

    obj = obj+(st-P(3:5))'*Q*(st-P(3:5))+delta_u'*Ru*delta_u +cons'*Qu*cons; % calculate obj

    st_next = X(:,k+1);

    

    valve =50^((con)-1);

  Pa = ((Ta*R/(Va*Mw)+(g*La)/Va).*x1);

  rho_a = ((Mw/(Ta*R))*Pa);

  Pwh = ((Tw*R/Mw)*(x2./(Lw*Aw-(x3./rho_o))));

  rho_w  = ((x2+x2)./(Lw*Aw));

  wpc = Cpc*sqrt(rho_w.*(max(0.0,(Pwh - Ps))))*valve;

  Pw = Pwh + ((g/Aw)*(x2+x3));

  Pbh = Pw + (rho_o*g*Lbh);

  wiv = Civ*sqrt(rho_a*(max (0.0,(Pa-Pw))));

  wpg = ((x2./max(0.0,(x2+x3))).*wpc);

  wpo = ((x3./max(0.0,(x2+x3)))*wpc);

  wro = Cr*sqrt(rho_o*(max(0.0,(Pr-Pbh))));

  wrg = (GOR*wro);

  

    

    f_value = f(st,con);

    st_next_euler = st+ (T*f_value);

    x1=st_next_euler(1);

    x2=st_next_euler(2);

    x3=st_next_euler(3);

    

    g1 = [g1;st_next-(st_next_euler)*T] ;% compute constraints

    g2 = [g2;delta_u];

%     g3 = [g3;rhs1(1,1)-wiv];

%     g4 = [g4;rhs1(2,1)-wro];

%     g5 = [g5;rhs1(3,1)-wrg];

%     g6 = [g6;rhs1(4,1)-wpg];

%     g7 = [g7;rhs1(5,1)-wpo];

        u_km1 = con;

%         x1-a

    a=x1+(wgl - wiv)*T

    b=x2+(wiv - wpg+ wrg)*T;

    c=x3+(wro - wpo)*T;

    

    x1=a;

    x2=b;

    x3=c;

    

end

g = [g1;g2];

% g = [g1;g2;g3;g4;g5;g6;g7];

% make the decision variable one column  vector

OPT_variables = [reshape(X,3*(N+1),1);reshape(U,N,1)];

 

nlp_prob = struct('f', obj, 'x', OPT_variables, 'g', g, 'p', P);

 

opts = struct;

opts.ipopt.max_iter = 2000;

opts.ipopt.print_level =0;%0,3

opts.print_time = 0;

opts.ipopt.acceptable_tol =1e-8;

opts.ipopt.acceptable_obj_change_tol = 1e-6;

 

% solver = nlpsol('solver', 'ipopt', nlp_prob,opts);

solver = nlpsol('solver', 'ipopt', nlp_prob,opts);

 

 

args = struct;

 

args.lbg1(1:3*(N+1)) = 0;

args.ubg1(1:3*(N+1)) = 0;

args.lbg2(1:(N)) = -du_max;

args.ubg2(1:(N)) = du_max;

% args.lbg3(1:(N+1)) = 0;

% args.ubg3(1:(N+1)) = 0;

% args.lbg4(1:(N+1)) = 0;

% args.ubg4(1:(N+1)) = 0;

% args.lbg5(1:(N+1)) = 0;

% args.ubg5(1:(N+1)) = 0;

% args.lbg6(1:(N+1)) = 0;

% args.ubg6(1:(N+1)) = 0;

% args.lbg7(1:(N+1)) = 0;

% args.ubg7(1:(N+1)) = 0;

args.lbg=[args.lbg1,args.lbg2];

args.ubg=[args.ubg1,args.ubg2];

 

% args.lbg=[args.lbg1,args.lbg2,args.lbg3,args.lbg4,args.lbg5,args.lbg6,args.lbg7];

% args.ubg=[args.ubg1,args.ubg2,args.ubg3,args.ubg4,args.ubg5,args.ubg6,args.ubg7];

 

args.lbx(1:3:3*(N+1),1) = x1min; %mga lower bound

args.ubx(1:3:3*(N+1),1) = x1max; %mga upper bound

args.lbx(2:3:3*(N+1),1) = x2min; %mgt lower bound

args.ubx(2:3:3*(N+1),1) = x2max; %mgt upper bound

args.lbx(3:3:3*(N+1),1) = x3min; % mot lower bound

args.ubx(3:3:3*(N+1),1) = x3max; %mot upper bound

 

args.lbx(3*(N+1)+1:3*(N+1)+N,1) = u_min; %u lower bound

args.ubx(3*(N+1)+1:3*(N+1)+N,1) = u_max; %u upper bound

 

 

t0 = 0;

x0 = [x10 ;x20; x30];    % initial condition.

xs = [1600 ; 7000; 200]; % Reference concentration.

u_km1 = 0.65;

udes=0.65;

xx(:,1) = x0; % xx contains the history of states

t(1) = t0;

 

u0 = 0.831*ones(N,1);        % control inputs 

X0 = repmat(x0,N+1,1)'; % initialization of the states decision variables

 

sim_tim =150; % total sampling times

 

% Start MPC*

mpciter = 0; 

xx1 = [];

u_cl=[];

cost_fun=[];

x_sp = [];

 

 

tic;

for k=1:sim_tim;

% while(norm((x0-xs),2) > 0.05 && mpciter < sim_tim / T)

    args.p   = [x0;xs;u_km1;udes]; % set the values of the parameters vector

       args.x0  = [reshape(X0,3*(N+1),1);u0]; % initial value of the optimization variables

    sol = solver('x0', args.x0, 'lbx', args.lbx, 'ubx', args.ubx,...

        'lbg', args.lbg, 'ubg', args.ubg,'p',args.p);

    u = reshape(full(sol.x(3*(N+1)+1:end)),1,N)'; % get controls only from the solution

    xx1(:,1:3,mpciter+1)= reshape(full(sol.x(1:3*(N+1)))',3,N+1)'; % get solution TRAJECTORY

    u_cl= [u_cl ; u(1,:)];

    t(mpciter+1) = t0;

    cost_fun=[cost_fun;full(sol.f)];

    % Apply the control and shift the solution

    [t0, x0, u0] = shift(T, t0, x0, u,f);

    x_sp = [x_sp, xs];

    xx(:,mpciter+2) = x0;

    X0 = reshape(full(sol.x(1:3*(N+1)))',3,N+1)'; % get solution TRAJECTORY

    % Shift trajectory to initialize the next step

    X0 = [X0(2:end,:);X0(end,:)];

    mpciter;

    mpciter = mpciter + 1;

     if mpciter>=75

        xs = [1800 ; 300; 8000];

        udes=0.95;

     end

     

    u_km1=u(1,:);

    

    x1=x1+(wgl - wiv)*T;

    x2=x2+(wiv - wpg+ wrg)*T;

    x3=x2+(wiv - wpg+ wrg)*T;

    

       end;

[OJ adukwu]

Good day house,

This is for discussion and not an announcement please. Show I repost or what do I do next?

Thanks.

OJ

[Bruno M L]

Hi!

You define f = Function('f',{states,controls},{rhs}); which means that f is a function of states (x1;x2;x3) and control (u). But there is no x or u in rhs.

This is one of the possible errors.

[ZHENG Zhuang]

Hi, Have you solved this problem? I also get the same error.

[Edson Sabino da Silva]

Hi, I have the same problem. My "rhs" equation contains all of the input variables and I am still not capable of creating the solver. 
Any help would be appreciated. 

[Edson Sabino da Silva]

Hi, everyone.

After reading similar posts in this group, I understood that I must tell the solver which inputs are variables or parameters.
To do so, it is necessary to make the parameters appear in the "P" vector and the variables in "OPT_variables" vector.

By doing this, the " are free" error is gone and I can now vary the parameters as the MPC loop goes, wich I was unable before knowing the function of the vector "P" in the solver creation/use