Re: Constructing optimizer

[Ali_E]

eh... I get this:
The number of cell elements in input does not match OPTIMIZER declaration
how to debug?

[Ali_E]

my code is:

clc;
clear;
close all;
%% Time structure
t0 = 0;
tf = 30;
dt = 1;
t = t0:dt:tf-dt;
%% System structure ===> x(k+1)=A*x(k)+B*u(k)
A = [1.02,-0.1;0.1,0.98];
B = [3.5;2.6];
G = [0.3,0;0,0.3];
C = [1,0];
%% MPC parameter
Q = 1;
R = 1;
N = 16;
M = blkdiag(Q,R);
w1 = randn(1,500);
w1 = w1-mean(w1);
w1 = w1/1;
muw1 = mean(w1);
sigmaw1 = var(w1);
w2 = randn(1,500);
w2 = w2-mean(w2);
w2 = w2/1;
muw2 = mean(w2);
sigmaw2 = var(w2);
w=[w1;w2];
sigmaw = [sigmaw1,0;0,sigmaw2];
v = randn(1,500);
v = v-mean(v);
v = v/1;
muv = mean(v);
sigmav = var(v);
x = randn(2,1);
y = C*x+v(:,1);
%mux{1} = [-5;2];
%muy{1} = C*mux{1};
h1y = [1;1];
h2u = 1;
h3u = -1;
epsilon = 0.01;
g1 = 0.1;
g2 = 1;
g3 = 1;
%% Pre-allocation
U1 = zeros(1,200);
mux1 = zeros(1,200);
mux2 = zeros(1,200);
f1 = cell(1,200);
Y1 = cell(1,200);
K1 = cell(1,200);
u2 = cell(1,200);
muy1 = zeros(1,200);
y1 = zeros(1,200);
ud1 = zeros(1,200);
x1 = zeros(1,200);
x2 = zeros(1,200);
%% Main loop
c=0;
for t=t0:dt:tf-dt
    c=c+1;
    F = sdpvar(repmat(1,1,N+1),repmat(1,1,N+1));    
    p = sdpvar(repmat(2,1,N-1),repmat(2,1,N-1));
    pN = sdpvar(1,1);
    muu = sdpvar(repmat(1,1,N),repmat(1,1,N)) ;
    mux = sdpvar(repmat(2,1,N),repmat(1,1,N));
    muy  = sdpvar(repmat(1,1,N),repmat(1,1,N));
    sigmay = sdpvar(repmat(1,1,N+1),repmat(1,1,N+1));
    sigmau = sdpvar(repmat(1,1,N),repmat(1,1,N));
    sigmax= sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));    
    constraint = [];
    objective = 0;
    for k=1:N-1
        mux{k+1} = A*mux{k}+B*muu{k};
        muy{k+1} = C*mux{k+1};
        objective = objective+1/2*(trace(M*p{k}));
        constraint1 = [sigmax{k+1},A*sigmax{k},B*F{k},G*sigmaw;(A*sigmax{k})',sigmax{k},zeros(2,1),zeros(2);(B*F{k})',zeros(1,2),sigmay{k},zeros(1,2);G*sigmaw,zeros(2,2),zeros(2,1),sigmaw]>=0;
        constraint2 = [sigmay{k+1},C*A*sigmax{k},C*B*F{k},C*G*sigmaw,sigmav;(C*A*sigmax{k})',sigmax{k},zeros(2,1),zeros(2,2),zeros(2,1);(C*B*F{k})',zeros(1,2),sigmay{k},zeros(1,2),zeros(1);(C*G*sigmaw)',zeros(2,2),zeros(2,1),sigmaw,zeros(2,1);sigmav,zeros(1,2),zeros(1),zeros(1,2),sigmav]>=0;
        constraint3 = [sigmau{k},(F{k});(F{k})',sigmay{k}]>=0;
        constraint4 = h1y'*mux{k}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1y'*(sigmay{k})*h1y);
        constraint5 = h2u'*muu{k}<=(1-(0.5*epsilon))*g2-((0.95/2*epsilon*g2*(1-0.95))*h2u'*sigmau{k}*h2u);
        constraint6 = h3u'*muu{k}<=(1-(0.5*epsilon))*g3-((0.95/2*epsilon*g3*(1-0.95))*h3u'*sigmau{k}*h3u);
        constraint7 = [p{k},[sigmay{k};F{k}],[muy{k};muu{k}];([sigmay{k};F{k}])',sigmay{k},zeros(1);([muy{k};muu{k}])',zeros(1),1]>=0;
        constraint = [constraint,constraint1,constraint2,constraint3,constraint4,constraint5,constraint6,constraint7];
    end
    objective=objective+(1/2*trace((trace(Q*pN))));
    constraint8 = h1y'*(mux{N})<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(0.1))*h1y'*(C*sigmax{N}*C'+sigmav)*h1y);
    constraint9 = [pN-(sigmay{k}),(muy{N});(muy{N})',1] >=0;
    constraint = [constraint,constraint8,constraint9];
    ops = sdpsettings('solver','mosek');
    %a = optimize(constraint,objective,ops);
    controller = optimizer(constraint,objective,ops,mux{1},muu{1});
    mux{1} = [-5;2];
    muy{1} = C*mux{1};
    %u = value(muu{1});
    u = controller{mux};
    mux3 = value(mux{1});
    mux1(c) = mux3(1,1);
    mux2(c) = mux3(2,1);
    muy3 = muy{1};
    F = value(F{1});
    Y = value(sigmay{1});
    x = value(x);
    y = value(y);
    U1(c) = u(1,1);
    f1{c} = F;
    Y1{c} = Y;
    K = F/Y;
    K1{c} = K;
    K2 = value(K1{c});
    u1 = u+K2*((y)-(muy3));
    u1 = value(u1);
    u2{c} = u1;
    x = A*x+B*u1+G*w(:,c);
    y = C*x+v(:,c);
    mux{1} = value(mux{2});
    sigmax{1} = value(sigmax{2});
    muy{1} = C*mux{1};
    sigmay{1} = C*sigmax{1}*C'+sigmav;
    muy1(c) = muy3;
    y1(c) = y;
    ud = value(u2{c});
    ud1(c) = ud(1,1);
    x1(c) = x(1,1);
    x2(c) = x(2,1);
end
%% Plot results
figure(1)
subplot(3,1,1)
plot(muy1,'linewidth',2)
xlabel('Time')
ylabel('muy1')
subplot(3,1,2)
plot(mux1,'linewidth',2)
xlabel('Time')
ylabel('mux1')
subplot(3,1,3);
plot(mux2,'linewidth',2)
xlabel('Time')
ylabel('mux2')

[Johan Löfberg]

% At this point mux{1} is an sdpvar vector
controller = optimizer(constraint,objective,ops,mux{1},muu{1});

mux{1} = [-5;2];
muy{1} = C*mux{1};
%u = value(muu{1

});
% Now you are sending a cell
u = controller{mux};

so makes no sense in terms of dimensions

Of course the code in general makes no sense as you are defining optimizer objects repeatedly thus only making the code slower compared to doing a call to optimize instead

[Ali_E]

ok according to https://yalmip.github.io/example/standardmpc/ it has to be something like:

clc;
clear;
close all;
%% Time structure
t0 = 0;
tf = 30;
dt = 1;

loop
F = sdpvar(repmat(1,1,N+1),repmat(1,1,N+1));

p = sdpvar(repmat(2,1,N-1),repmat(2,1,N-1));
pN = sdpvar(1,1);
muu = sdpvar(repmat(1,1,N),repmat(1,1,N)) ;
mux = sdpvar(repmat(2,1,N),repmat(1,1,N));
muy  = sdpvar(repmat(1,1,N),repmat(1,1,N));
sigmay = sdpvar(repmat(1,1,N+1),repmat(1,1,N+1));
sigmau = sdpvar(repmat(1,1,N),repmat(1,1,N));
sigmax= sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));

constraints = [];

objective = 0;
for k=1:N-1
    mux{k+1} = A*mux{k}+B*muu{k};
    muy{k+1} = C*mux{k+1};
    objective = objective+1/2*(trace(M*p{k}));
    constraint1 = [sigmax{k+1},A*sigmax{k},B*F{k},G*sigmaw;(A*sigmax{k})',sigmax{k},zeros(2,1),zeros(2);(B*F{k})',zeros(1,2),sigmay{k},zeros(1,2);G*sigmaw,zeros(2,2),zeros(2,1),sigmaw]>=0;
    constraint2 = [sigmay{k+1},C*A*sigmax{k},C*B*F{k},C*G*sigmaw,sigmav;(C*A*sigmax{k})',sigmax{k},zeros(2,1),zeros(2,2),zeros(2,1);(C*B*F{k})',zeros(1,2),sigmay{k},zeros(1,2),zeros(1);(C*G*sigmaw)',zeros(2,2),zeros(2,1),sigmaw,zeros(2,1);sigmav,zeros(1,2),zeros(1),zeros(1,2),sigmav]>=0;
    constraint3 = [sigmau{k},(F{k});(F{k})',sigmay{k}]>=0;
    constraint4 = h1y'*mux{k}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1y'*(sigmay{k})*h1y);
    constraint5 = h2u'*muu{k}<=(1-(0.5*epsilon))*g2-((0.95/2*epsilon*g2*(1-0.95))*h2u'*sigmau{k}*h2u);
    constraint6 = h3u'*muu{k}<=(1-(0.5*epsilon))*g3-((0.95/2*epsilon*g3*(1-0.95))*h3u'*sigmau{k}*h3u);
    constraint7 = [p{k},[sigmay{k};F{k}],[muy{k};muu{k}];([sigmay{k};F{k}])',sigmay{k},zeros(1);([muy{k};muu{k}])',zeros(1),1]>=0;

    constraints = [constraints,constraint1,constraint2,constraint3,constraint4,constraint5,constraint6,constraint7];

end
objective=objective+(1/2*trace((trace(Q*pN))));
constraint8 = h1y'*(mux{N})<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(0.1))*h1y'*(C*sigmax{N}*C'+sigmav)*h1y);
constraint9 = [pN-(sigmay{k}),(muy{N});(muy{N})',1] >=0;

constraints = [constraints,constraint8,constraint9];

ops = sdpsettings('solver','mosek');
%a = optimize(constraint,objective,ops);

controller = optimizer(constraints,objective,ops,mux{1},muu{1});
%mux = {[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2],[-5;2]};

mux{1} = [-5;2];
muy{1} = C*mux{1};
%u = value(muu{1});

u = controller{mux{1}};

c=0;
for t=t0:dt:tf-dt
    c=c+1;

    mux3 = value(mux{1});
    mux1(c) = mux3(1,1);
    mux2(c) = mux3(2,1);
    muy3 = muy{1};
    F = value(F{1});
    Y = value(sigmay{1});
    x = value(x);
    y = value(y);
    U1(c) = u(1,1);
    f1{c} = F;
    Y1{c} = Y;
    K = F/Y;
    K1{c} = K;
    K2 = value(K1{c});
    u1 = u+K2*((y)-(muy3));
    u1 = value(u1);
    u2{c} = u1;
    x = A*x+B*u1+G*w(:,c);
    y = C*x+v(:,c);
    mux{1} = value(mux{2});
    sigmax{1} = value(sigmax{2});

    muy{1} = C*mux{1};

    sigmay{1} = C*sigmax{1}*C'+sigmav;
    muy1(c) = muy3;
    y1(c) = y;
    ud = value(u2{c});
    ud1(c) = ud(1,1);
    x1(c) = x(1,1);
    x2(c) = x(2,1);
end
%% Plot results
figure(1)
subplot(3,1,1)
plot(muy1,'linewidth',2)
xlabel('Time')
ylabel('muy1')
subplot(3,1,2)
plot(mux1,'linewidth',2)
xlabel('Time')
ylabel('mux1')
subplot(3,1,3);
plot(mux2,'linewidth',2)
xlabel('Time')
ylabel('mux2')

but again it gives unfamiliar error:
Brace indexing is not supported for variables of this type.

Error in Optimizer_New (line 104)

[Johan Löfberg]

to begin with, value operator can never be used when using optimizer objects. all solutions are communicated via the outputs from the optimizer call , and secondly you never use your optimizer object in the loop, meaning you really don't understand the concept of it yet. read the wiki example more clearly to see how it is used

[Ali_E]

value operator can never be used when using optimizer objects ===> so how we extract results and address them and define them and plot them etc?
you never use your optimizer object in the loop ===> what do you mean "optimizer object"? you mean "mux{1}"?

[Johan Löfberg]

P = optimizer(constraints,objective,options,{parameter1,parameter2,...},{output1,output2,...})

...

parametervalue1 = someconstant
...

for ...

  solution = P(parametervalue1,parametervalue2,...)

  theoptimalvalueofoutput1 = solution{1};
  theoptimalvalueofoutput2 = solution{2};

  parametervalue1 = ...

end
...

[Ali_E]

question here is that how should I address value(muu{1}) when I use optimizer because when I use optimize value of muu{1} is a cell and I address to that, but when I use optimizer I donnow how to address that muu{1} because I defined output as muu{:}, but I don't want to put muu{:} in close loop only muu{1} well be used and brace indexing is not allowed.

[Ali_E]

P = optimizer(constraints,objective,options,{parameter1,parameter2,...},{output1,output2,...})

...

parametervalue1 = someconstant
...



what are these parameter1,parameter2,... and how should I choose them???
and parametervalue1 is initial value for parameter1?

[Johan Löfberg]

You don't use value, but declare muu{1} as one of the output parameters when you create the optimizer object. muu{1} is not a cell, it is a scalar in your code.

[Johan Löfberg]

I don't know what the parameters are in your model, because you've still haven't managed to explain what is changing in the iterations, i.e. what the cause for solving new optimization problems is because sometimes when you show your code, there is nothing which changes so the use of an optimzer approach is moot). In almost all MPC scenarios, there is only one parametric variable (the thing that the optimization problems is parameterized in), and that is the initial state x{1} in the notation of the MPC example, i.e. the MPC problem you solve looks completely the same the whole time, except that the initial state changes. 

[Ali_E]

I need muy and muu and muy{1} = C*mux{1} and I want to give them to an affine feedback u1 = u+K2*((y)-(muy3)) and get u1 to give to x = A*x+B*u1+G*w(:,c) and get x to give to y = C*x+v(:,c) and get y

[Johan Löfberg]

trying to decipher that incoherent sentence and assuming muy is a decision variable and x is a parameter and muy3 means muy{1}(3) , you are saying

P = optimizer(constraints,objective,options,x{1},{u,muy{1}})

x = initialvalue

for loop

 sol = P(x);

 y = C*x+v(:,c);

 q = C*sol{2};; 

 u1 = sol{1}+K2*((y)-q(3))

 x = A*x + B*u1 + G*w(:,c);

end

[Ali_E]

excuse me professor it's not that easy I got a bunch of constraints and there are relationship between some parameters in those constraints.

for t=t0:dt:tf-dt
    c=c+1;

    sigmax= sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));

    sigmay = sdpvar(repmat(1,1,N+1),repmat(1,1,N+1));
    sigmau = sdpvar(repmat(1,1,N),repmat(1,1,N));

    p = sdpvar(repmat(2,1,N-1),repmat(2,1,N-1));
    pN = sdpvar(1,1);

    F = sdpvar(repmat(1,1,N+1),repmat(1,1,N+1));
    %mux = sdpvar(repmat(2,1,N),repmat(1,1,N));
    muu = sdpvar(repmat(1,1,N),repmat(1,1,N)) ;
    %muy  = sdpvar(repmat(1,1,N),repmat(1,1,N));
    constraint = [];

    objective = 0;
    for k=1:N-1
        mux{k+1} = A*mux{k}+B*muu{k};
        muy{k+1} = C*mux{k+1};
        objective = objective+1/2*(trace(M*p{k}));

        constraint2 = [sigmay{k+1},C*A*sigmax{k},C*B*F{k},C*G*sigmaw,sigmav;(C*A*sigmax{k})',sigmax{k},zeros(2,1),zeros(2,2),zeros(2,1);(C*B*F{k})',zeros(1,2),sigmay{k},zeros(1,2),zeros(1);(C*G*sigmaw)',zeros(2,2),zeros(2,1),sigmaw,zeros(2,1);sigmav,zeros(1,2),zeros(1),zeros(1,2),sigmav]>=0;
        constraint1 = [sigmax{k+1},A*sigmax{k},B*F{k},G*sigmaw;(A*sigmax{k})',sigmax{k},zeros(2,1),zeros(2);(B*F{k})',zeros(1,2),sigmay{k},zeros(1,2);G*sigmaw,zeros(2,2),zeros(2,1),sigmaw]>=0;
        constraint3 = [sigmau{k},(F{k});(F{k})',sigmay{k}]>=0;
        constraint4 = h1y'*mux{k}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1y'*(sigmay{k})*h1y);
        constraint5 = h2u'*muu{k}<=(1-(0.5*epsilon))*g2-((0.95/2*epsilon*g2*(1-0.95))*h2u'*sigmau{k}*h2u);
        constraint6 = h3u'*muu{k}<=(1-(0.5*epsilon))*g3-((0.95/2*epsilon*g3*(1-0.95))*h3u'*sigmau{k}*h3u);
        constraint7 = [p{k},[sigmay{k};F{k}],[muy{k};muu{k}];([sigmay{k};F{k}])',sigmay{k},zeros(1);([muy{k};muu{k}])',zeros(1),1]>=0;

        constraint = [constraint,constraint1,constraint2,constraint3,constraint4,constraint5,constraint6,constraint7];

    end
    objective=objective+(1/2*trace((trace(Q*pN))));
    constraint8 = h1y'*(mux{N})<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(0.1))*h1y'*(C*sigmax{N}*C'+sigmav)*h1y);
    constraint9 = [pN-(sigmay{k}),(muy{N});(muy{N})',1] >=0;

    constraint = [constraint,constraint8,constraint9];

and also that K2 in u1 = u+K2*((y)-(muy3)) which evaluates from (using optimize):

F = value(F{1});
    Y = value(sigmay{1});

    f1{c} = F;
    Y1{c} = Y;
    K = F/Y;
    K1{c} = K;
    K2 = value(K1{c});

so F and sigmay are outputs too, but the main question is that only decision variables are parameters in P = optimizer(constraints,objective,options,{parameter1,parameter2,...},{output1,output2,...})

and another question is that initial constant values are arbitrary or not???


[Johan Löfberg]

Although that is entirely inconsistent with the code you have. In the code, it appears mux{1} is the parameter having some kind of role of the state (x is never used) and it also looks like both F and Y are needed to compute the input as they define K, which is new in every iteration. Hence more likely something along the lines of 

P= optimizer(constraint,objective,options,mux{1},{muu{1},F{1},sigmay{1}})...

[Ali_E]

yeah it might be something like this, but there are two questions:
1. what is parameter's definition? is it decision variable or not? what is the exact meaning of parameter in:

P = optimizer(constraints,objective,options,{parameter1,parameter2,...},{output1,output2,...})

2. are initial values (that I should declare after P=optimizer(..)) arbitrary or not?

[Johan Löfberg]

I still doubt you've read these. You are solving exactly the same problem as illustrated in the MPC example, except that you have several outputs from the optimizer object to compute u.

https://yalmip.github.io/example/standardmpc/

https://yalmip.github.io/command/optimizer/

[Ali_E]

eh... I've read them, but there is no definition for "parametee1,parameter2,..." there
and also there is nothing about initial values and also any answer to "are initial values arbitrary or not?"

[Johan Löfberg]

Of course there are

In the first example, the current state is the only parameter, and the whole control sequence is the output

controller = optimizer(constraints, objective,[],x{1},[u{:}]);

The second example involves a more advanced parameterization of an optimization problem (current state, preview of reference, estimated disturbance, previous input,  and B matrix in the simulation model) and outputs whole input and state trajectory

parameters_in = {x{1},[r{:}],d,pastu,B};
solutions_out = {[u{:}], [x{:}]};

controller = optimizer(constraints, objective,sdpsettings('solver','gurobi'),parameters_in,solutions_out);

Initial values: Well that is your simulation. You are simulating a dynamical system from some initial state.and it is up to you to select that initial state

[Ali_E]

ok I tried this, but still an error shows:

parameters_in = {F{1},p{1},pN,muu{1},mux{1},muy{1},sigmay{1},sigmau{1},sigmax{1}};
solutions_out = {[mux{1}],[muu{1}],[F{1}],[sigmay{1}],[muy{1}]};
controller = optimizer(constraints,objective,ops,parameters_in,solutions_out);
F{1} = 0;
muu{1} = 0;

mux{1} = [-5;2];

muy{1} = C*mux{1};
sigmax{1} = [1,0;0,1];

sigmay{1} = C*sigmax{1}*C

'+sigmav;

c=0;
for t=t0:dt:tf-dt
    c=c+1;

    sol = controller(F{1},p{1},pN,muu{1},mux{1},muy{1},sigmay{1},sigmau{1},sigmax{1});
    mux3 = sol{1};
    u = sol{2};
    F = sol{3};
    Y = sol{4};
    muy3 = sol{5};

    mux1(c) = mux3(1,1);
    mux2(c) = mux3(2,1);

    x = value(x);
    y = value(y);
    U1(c) = u(1,1);
    f1{c} = F;
    Y1{c} = Y;
    K = F/Y;
    K1{c} = K;
    K2 = value(K1{c});
    u1 = u+K2*((y)-(muy3));
    u1 = value(u1);
    u2{c} = u1;
    x = A*x+B*u1+G*w(:,c);
    y = C*x+v(:,c);

    sol{1} = value(mux{2});

    sigmax{1} = value(sigmax{2});

    sol{5} = C*sol{1};

    sigmay{1} = C*sigmax{1}*C'

+sigmav;

    muy1(c) = muy3;
    y1(c) = y;
    ud = value(u2{c});
    ud1(c) = ud(1,1);
    x1(c) = x(1,1);
    x2(c) = x(2,1);
end
%% Plot results
figure(1)
subplot(3,1,1)
plot(muy1,'linewidth',2)
xlabel('Time')
ylabel('muy1')
subplot(3,1,2)
plot(mux1,'linewidth',2)
xlabel('Time')
ylabel('mux1')
subplot(3,1,3);
plot(mux2,'linewidth',2)
xlabel('Time')
ylabel('mux2')

error: Undefined function 'isnan' for input arguments of type 'sdpvar'.

[Johan Löfberg]

code makes absolutely no sense. 

why is F{1} a parameter and an output. same with muu and sigmay

why are you bracketing like this [mux{1}]

 x = value(x); huh? x is a numerical value so value makes no sense

. same with K2 = value(K1{c});

I really really think your only parameter is mux{1}, which I guess is the name for the expected value of state x (hence the expected value of current state is of course the given value of x, so you would be making the call controller(x)

[Johan Löfberg]

and here is my guess what you are doing

You are parameterizing the input along the trajectory as something like muu{k} + K(y(k) - Cmux(k)), i.e. a fixed value +  feedback of deviation between actual output y(k) and expected output.

However, this is MPC, so the problem is repeatedly resolved, so the solutions along the trajectory are never used, they are just used to estimate the optimal worst-case cost. Hence, the only actual input will be muu{1} (since there will be no feedback active at the initial state, as the current state is known and thus K(y(k) - Cmux(k)) = K(y(k) - Cx(k))  = 0

Hence, optimizer has parameter mux{1} and output muu{1}, and the code is

controller = optimizer(constraints,objective,ops,mux{1},muu{1}

);
x =
for t = ...
 u =controller(x);
 x =A*x+B*u+...
end

i.e., standard MPC simulation

[Ali_E]

let me define it line by line:
I need these values as solutions:

solutions_out = {[mux{1}],[muu{1}],[F{1}],[sigmay{1}],[muy{1}]};

so I have:

mux3 = sol{1};
    u = sol{2};
    F = sol{3};
    Y = sol{4};
    muy3 = sol{5};

then I define my time loop:

for t=t0:dt:tf-dt
    c=c+1;

    %sol = controller(F{1},p{1},pN,muu{1},mux{1},muy{1},sigmay{1},sigmau{1},sigmax{1});
    sol=controller(mux{1});

    mux3 = sol{1};
    u = sol{2};
    F = sol{3};
    Y = sol{4};
    muy3 = sol{5};
    mux1(c) = mux3(1,1);
    mux2(c) = mux3(2,1);

    U1(c) = u(1,1);
    f1{c} = F;
    Y1{c} = Y;
    K = F/Y;
    K1{c} = K;
    K2 = value(K1{c});
    u1 = u+K2*((y)-(muy3));
    u1 = value(u1);
    u2{c} = u1;
    x = A*x+B*u1+G*w(:,c);
    y = C*x+v(:,c);
    sol{1} = value(mux{2});
    sigmax{1} = value(sigmax{2});
    sol{5} = C*sol{1};
    sigmay{1} = C*sigmax{1}*C

'+sigmav;
    muy1(c) = muy3;
    sol{4} = C*sigmax{1}*C'+sigmav;

    y1(c) = y;
    ud = value(u2{c});
    ud1(c) = ud(1,1);
    x1(c) = x(1,1);
    x2(c) = x(2,1);
end

muu and F and sigmax and p are my decision variables. mux acts as state. there are also other parameters in constraints so I have to define them as cells using repmat:

%parameters_in = {F{1},p{1},pN,muu{1},mux{1},muy{1},sigmay{1},sigmau{1},sigmax{1}};

solutions_out = {[mux{1}],[muu{1}],[F{1}],[sigmay{1}],[muy{1}]};

controller = optimizer(constraints,objective,ops,mux{1},solutions_out);

I also attached the file that I want to construct optimizer for it.

[Johan Löfberg]

If mux is supposed to represent the expected value of the states, why are you initializing it to [-5; 2] and the true state to randn, i.e. absolutely no connection

[Ali_E]

muy=C*mux then u1 = u+K2*((y)-(muy3))  and muy3=muy{1} and u=muu{1} and at last x = A*x+B*u1+G*w(:,c) so there is connection, but that is not straight connection etc

[Johan Löfberg]

You initialize them to something completely arbitrary vs a very specific value. That does not indicate that there is any connection.

Does mux represent the expected value of x when you move along the trajectory disturbed by external disturbance. Then you would have mux{1} = x

Alternatively, if mux{1} represents a state estimate of the true state x, and that state estimate and the control sequence, is computed by the optimization problem (thus mux{2} being the new estimate), it would make sense to initialize x to some value, and mux{1} to x + random initial estimation error

[Ali_E]

mux is the mean of state x

[Johan Löfberg]

yes, but is x measured (and thus mux{1}=x) or are you assuming that you use mux as the estimator (which would make more sense here as you are parameterizing the input in terms of measurements. If you had the state available for the controller, you would parameterize in x)

[Ali_E]

x is not measured. whole problem is this that x is not measured. x is stochastic and thus we have an stochastic problem. we convert it to a deterministic problem. now we have mux and a deterministic problem. this was the whole idea.

[Johan Löfberg]

and thus the question is if you use the solution of mux{2} as your state estimator update

[Ali_E]

yeah there is sol{1} = value(mux{2}) and sol{1}=value(mux{1}) in my code

[Johan Löfberg]

and thus mux{1} would be the parameter, and mux{2}, muu{1}, F{1}, sigmay{1} and muy{1} would be the outputs

...
...
Controller = optimizer(con,obj,ops,mux{1},{mux{2},muu{1},F{1}, sigmay{1} and muy{1} )
x = randn(2,1)
mux1 = [-2;5]% or what ever you logic is for initial guess of estimate ...
for t = ...
sol = Controller(mux1),
construct u1 from solution
x = A*x + B*u1 + ...
% Take mux{2} solution as new state estimate
mux1 = sol{1};
end