How to model constraint on the shape of the decision variables?

[Hongkai Chen]

Hi Johan,

I am modeling a state estimator using moving horizon estimation. Instead of making all disturbances decision variables in the optimization problem, I have the information that in the window I will have three square-wave disturbances.

I want to use three amplitudes, durations and start times to replace window size number of decision variables. I also have the constraints on the amp, duration and time too.

I know in the MHE I have to make the disturbances a window size long sdpvar, and use it in the for-loop of the description of the dynamic system. But I don't know how to do it if I want to constrain on the shape and number of disturbances.

Thank you!

[Johan Löfberg]

sounds like a combinatorial constraint, i.e. you need binary variables and logic constraints for this, similiar to minimum up-time etc

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

[Hongkai Chen]

For example, if the code is like this

A = rand(6,6);

b = rand(6,1);

C = rand(6,1):

D =rand(6,1);

E = rand(6,1);

u = sdpvar(30,1);

w = sdpvar(30,1);

x = sdpvar(6,1);

x_hat = sdpvar(6,1);

y= sdpvar(30,1);

Y = [];

V=[];

xk = x;

for k = 1:30

    

    xk =  A*xk + b + C*u(k)  +D*w(k) ;

    Y = [Y; E*xk];

    V = [V; y(k) - E*xk];

end

F = [0 <= u <= 10]; 

objective = norm(x - x_hat,1) + norm(V,1);

G = [ 0 <= w <= 20];       

ops = sdpsettings('solver','mosek','verbose',0);

controller = optimizer(G,objective,ops,{u, x_hat, y},  {xk,w});

Now I have more information about w, e.g. I know there should be one or two square waves, with the duration between 3 and 5 and amplitude between 10 and 15. How do I model this?

[Johan Löfberg]

w looks to be an uncertainty, so to have it as a parameter in an optimizer object makes no sense. 

if it is an uncertainty, you have to define how you want to cope with it. Worst-case? In any way, I have no idea how o would robustify against a disturbance with your character, sounds extremely complicated

[Hongkai Chen]

Actually, this is a state estimation. so w becomes decision variables in this problem. Now I only have simple ub and lb on it.

[Johan Löfberg]

ok, then you should be able to use the minup/mindown logic models as a template

[Hongkai Chen]

Hi professor,

I used the model as a template as following:

A = rand(6,6);

b = rand(6,1);

C = rand(6,1);

D = rand(6,1);

N_SE = 150; % estiamtion window size

% dimension

nu_SE = 1;

nw_SE = 1; % dimension of disturabances

nx_SE = 6; % dimension of states

% The model will be composed of two variables. The first is a binary variable which 

% controls if a disturbance is taken at time k, and the second variable is a continuous 

% variable representing the disturbance amount taken at time k.

onoff = binvar(N_SE, nw_SE);

w_SE = sdpvar(N_SE,nw_SE); % Disturbance w(0),...,w(N_SE - 1)

% variables, now the decision variables are the states and the disturbances

x_SE = sdpvar(nx_SE,1); % Simulated initial States x(0), this is a decision variable

u_SE = sdpvar(nu_SE, N_SE); % insulin input u(0),...,u(N_SE - 1). here they are known values.

x_hat_SE = sdpvar(nx_SE,1); % Estimated initial states x_hat(0)

y_SE = sdpvar( 1, N_SE); % measurements y(0),...,y(N_SE - 1)

% system dynamics

E = [1 0 0 0 0 0];%times xk to get the first parameter

Y_SE = [];% Y: system output % simulated output H(x(t))

V_SE = [];% V: discrepency between measurements and simulated outputs

xk_SE = x_SE;

for k = 1:N_SE

    xk_SE = A*xk_SE + b + C*u_SE(k)  + D*w_SE(k);

    Y_SE = [Y_SE; E*xk_SE];

    V_SE = [V_SE; y_SE(k) - E*xk_SE];

end

% The amount of meals are either between the ub and lb (when onoff key is on)

% , or being 0 (when onoff key is off).

upperbound = sdpvar(N_SE,nw_SE);

lowerbound = sdpvar(N_SE,nw_SE);

Constraints =[];

for k = 1:N_SE

    Constraints = [Constraints, onoff(k)*lowerbound(k) <= w_SE(k)<= onoff(k)*upperbound(k)];

end

% Constraints that limit the number of disturbance in a window

limit = 3;

indicator = zeros(N_SE-1,1);

for k = 2:N_SE

  % indicator will be 1 only when switched on

  indicator(k-1) = onoff(k)-onoff(k-1);

end

% At most limit number of 1's

Constraints = [Constraints, nnz(indicator == 1) <= meal_limit, nnz(indicator == 1) == nnz(indicator == -1) ];

% disturbances are square, so all the values between a (1, -1) pair are the same

meal_number = nnz(indicator == 1);

[s, location] = sort(indicator);

if meal_number > 0

    for k = 1:meal_number

        range = location(N_SE - meal_number -1 + k):location(k)-1;

        Constraints = [Constraints, w_SE(range) == w_SE(range(1))];

    end

end

objective = norm(x_SE - x_hat_SE, 1) + norm(V_SE, 1); 

ops = sdpsettings('verbose',2,'debug',1);

% use optimizer instead of optimize

estimator = optimizer(Constraints, objective, ops, {u_SE, y_SE, x_hat_SE,upperbound,lowerbound}, {xk_SE, w_SE, onoff,objective} );

My assumption is in a moving horizon estimation window, there could be a limited number of square-wave disturbances. However, I keep getting 0.5 for all onoff values? What is the problem of my model? How can I fix it with the assumption? And what is the best solver for this?

Thank you very much!

[Johan Löfberg]

Undefined function or variable 'meal_limit'.

 

[Johan Löfberg]

and you never ever multply binary variables and other variables

onoff(k)*lowerbound(k)

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

[Johan Löfberg]

ah, bounds are parameters in optimizer, then it should work with the multiplication

[Johan Löfberg]

and your code is not particularly interesting to debug since you never use the optimizer object, hence we have no idea what happens when you use it. 

Nor do we know which solver you are using. I guess it will select bnb, as the problem is a nonconvex nonlinear program when the bounds are parameter and thus initially are decision variables when YALMIP analyzes the model. YALMIP has to be told explicitly which solver to use if you want to override its analysis, with you promising that the model actually is a MILP once the parameters are fixed.

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

[Johan Löfberg]

and this

% At most limit number of 1's

Constraints = [Constraints, nnz(indicator == 1) <= meal_limit, nnz(indicator == 1) == nnz(indicator == -1) ];

% disturbances are square, so all the values between a (1, -1) pair are the same

meal_number = nnz(indicator == 1);

[s, location] = sort(indicator);

if meal_number > 0

   for k = 1:meal_number

       range = location(N_SE - meal_number -1 + k):location(k)-1;

       Constraints = [Constraints, w_SE(range) == w_SE(range(1))];

   end

end

is going to be an absolutely massive MILP model. The sort operator alone will introduce over 20000 binary variables. The nnz operators are probably pretty expensive too. Surely you must be able to express the fact that the onoff variable is a pulse in a much much compact and cheaper way.

[Johan Löfberg]

the code shouldn't even be possible to run

meal_number= nnz(indicator == 1);

[s,location] = sort(indicator);

if meal_number > 0

meal_number is a decision variable, and those can never be used in if-statements

[Hongkai Chen]

Sorry, that is just 'limit'.

[Hongkai Chen]

what is the best solver to solve this problem?

[Hongkai Chen]

Could you give me any idea how can I express the pulses in a better way?

[Johan Löfberg]

it looks MILP representable, so gurobi/mosek/cplex 

[Johan Löfberg]

Your code is not valid as indicator will be a vector of nans due to your initialization as zeros

indicator = zeros(N_SE-1,1);

for k = 2:N_SE

  % indicator will be 1 only when switched on

  indicator(k-1) = onoff(k)-onoff(k-1);

end

with that, the nnz operators simply returns 0

https://yalmip.github.io/naninmodel/

[Hongkai Chen]

I told YALMIP to use gurobi. The results came out in a second. The onoff are binary now. But the disturbance are definitely not pulse. so I guess the express of pulses are fine, but the if-statement is the problem?

Academic license - for non-commercial use only

Optimize a model with 614 rows, 464 columns and 24480 nonzeros

Variable types: 314 continuous, 150 integer (150 binary)

Coefficient statistics:

  Matrix range     [5e-07, 2e+01]

  Objective range  [1e+00, 1e+00]

  Bounds range     [1e+00, 1e+00]

  RHS range        [2e-03, 5e+04]

Presolve removed 460 rows and 158 columns

Presolve time: 0.02s

Presolved: 154 rows, 306 columns, 11783 nonzeros

Variable types: 306 continuous, 0 integer (0 binary)

Root relaxation: objective 4.237454e+01, 440 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0      42.3745350   42.37454  0.00%     -    0s

Explored 0 nodes (440 simplex iterations) in 0.06 seconds

Thread count was 4 (of 4 available processors)

Solution count 1: 42.3745 

Optimal solution found (tolerance 1.00e-04)

Best objective 4.237453504107e+01, best bound 4.237453504107e+01, gap 0.0000%

xk_next =

   1.0e+04 *

    0.0008

    0.0192

    0.0219

    0.0060

    5.2745

    0.0009

estimated_disturbance =

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

    7.6539

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

   18.3607

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

   20.8333

    2.6178

   20.8333

   19.3379

   17.6999

   17.6950

   17.7046

   17.7094

   13.4749

    3.1115

    2.4347

    2.8359

    3.2628

    3.6291

    3.9433

    4.2129

    4.4454

    4.6607

    4.6833

    4.8525

    4.9962

    5.1181

    5.2215

    5.3092

    2.4569

   17.1341

         0

         0

         0

         0

         0

         0

         0

         0

         0

         0

    8.0971

   20.8333

    6.2142

    6.6353

    7.7290

    7.6640

    7.6054

    7.5531

    7.5067

    7.4659

    7.4307

    7.4047

    7.3405

    7.3248

    7.3114

    7.3002

    7.2911

    7.2840

    7.2788

    7.2753

    7.2742

    7.2829

    7.2144

    7.2291

    7.2423

    7.2544

    7.2657

    7.2764

    7.2868

    7.2970

    7.3073

    6.8590

   13.0720

         0

         0

         0

         0

         0

         0

    9.3122

    7.7629

    5.7254

    6.1284

    6.0784

    6.0350

    5.9973

    5.9645

    5.9360

    5.9112

    5.8896

    5.8708

    5.8544

    5.8391

    5.8266

    5.8156

    5.8058

    5.7971

    5.7893

    5.7823

    5.7760

         0

         0

onoff =

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

     1

objective_value =

   42.3745

[Johan Löfberg]

as I said, your indicator vector is nonsense

[Hongkai Chen]

So the solution is 

indicator = [];

for k = 2:N_SE

  % indicator will be 1 only when switched on

  indicator = [ indicator onoff(k)-onoff(k-1)];

end

Is this correct?

[Johan Löfberg]

Think in terms of if-statements, and write those as implies (or manually model logic as here to understand what you actually do)

Allow at most two pulses, of length 3 to 5

if x(i)=1 and x(i-1) = 0 then x(i:i+2)=1,x(i+5)=0  % Step, so at least three coming are 1, but no more than 5

if x(i)=1 and x(i-1)=0 and x(i+j) = 1 then x(i:j) = 1 % You cannot have a whole in a pulse, for j = 2:4

Convert

x(i:i+2) >= x(i) + (1-x(i-1))-1

1-x(i+5) >= x(i) + (1-x(i-1)) -1

x(i:j) >= x(i) + (1-x(i-1))+x(i+j) -2

or something like that to model a pulse

Nwo you have to say that you have at most 2 pulses. Either introduce indicator variables (which would be used above) and constrain the sum of those, or sum(max(diff(x),0)<= 2

and then you have to deal with the special case that you have a pulse that starts already t i = 1 (and thus no step up) and pulses going all the way to the end (and thus never move back to 0)

[Johan Löfberg]

something like that yes would be the brute version, the matlab standard would be indicator = diff(onoff)

[Johan Löfberg]

or better yet, reuse the code from the unit commitment example better

If you have at most three pulses, then write pulse = pulse1 + pulse2 + pulse3 where each vector is binary

Add constraints for minup(pulse1), minup(pulse2),minup(pulse3)

Pulses have to be separated so pulse1+pulse2+pulse3 <= 1

You also have upper bound on pulse length which simply is sum(pulse1) <= nupper, sum(pulse3) <= nupper,sum(pulse3) <= nupper

[Johan Löfberg]

something like this

n = 20;

y = rand(20,1) > .4;

x1 = binvar(n,1);

x2 = binvar(n,1);

x = x1 + x2;

Model = [x1 + x2 <= 1]

C1 = consequtiveON([x1],3);

Model = [Model, C1, sum(max(diff([0;x1]),0)) <= 1, sum(x1) <= 6];

C2 = consequtiveON([x2],3);

Model = [Model, C2, sum(max(diff([0;x2]),0)) <= 1, sum(x2) <= 6];

optimize(Model,norm(y-x,1));

value([y x1 x2 x])

with changes required if you want to treat the borders differently etc

[Hongkai Chen]

If I also want to say the pulses are square and within boundaries. Can I do this?

n = 20;

y = rand(20,1) > .4;

x1 = binvar(n,1);

x2 = binvar(n,1);

x = x1 + x2;

 

Model = [x1 + x2 <= 1]

C1 = consequtiveON([x1],3);

Model = [Model, C1, sum(max(diff([0;x1]),0)) <= 1, sum(x1) <= 6];

C2 = consequtiveON([x2],3);

Model = [Model, C2, sum(max(diff([0;x2]),0)) <= 1, sum(x2) <= 6];

ub = sdpvar(n,1);

lb = sdpvar(n,1);

sdpvar z;

for i = 1:n

  Model =[Model, x(i)*lb(i) <=z <=x(i)*ub(i)];

end

 

optimize(Model,norm(y-z,1));

[Johan Löfberg]

x simply defines the structure, and you can use that variable further to multiply it with what ever you want etc

[Hongkai Chen]

With the latest modification, the pulse contained in boundaries. The problem is that pulses have to be square. 

I guess I cannot do something like this.

sdpvar a,b

x = x1*a + x2*b;

Because we cannot multiple binvar with sdpvar. Is there a way to do this?

[Johan Löfberg]

you use precisely the same models as in the unit commitment examples to represent products

x = z1+z2

if x1 then z1=a

if ~x1 then z1=0

etc

and those you model using logics, i.e. either implies(x1,z1==a), or manually through big-M, -Mx1 <= z-1a <= M*x1 etc.

[Hongkai Chen]

A = rand(6,6);

b = rand(6,1);

C = rand(6,1);

D = rand(6,1);

% parameter definition

N_SE = 150; % estiamtion window size

% dimension

nu_SE = 1;

nw_SE = 1; % dimension of disturabances

nx_SE = 6; % dimension of states

Constraints =[];

distur_1 = binvar(N_SE, nw_SE);

distur_2 = binvar(N_SE, nw_SE);

distur_3 = binvar(N_SE, nw_SE);

w1 = sdpvar(N_SE, nw_SE);

w2 = sdpvar(N_SE, nw_SE);

w3 = sdpvar(N_SE, nw_SE);

sdpvar distur_1_amp distur_2_amp distur_3_amp;

w_SE = w1 + w2 + w3; % Disturbance w(0),...,w(N_SE - 1)

Constraints = [Constraints, implies(distur_1,w1==distur_1_amp),implies(distur_2,w2==distur_2_amp),implies(distur_3,w3==distur_3_amp)];

Constraints = [Constraints, implies(~distur_1,w1==0),implies(~distur_2,w2==0),implies(~distur_3,w3==0)];

% variables, now the decision variables are the states and the disturbances

x_SE = sdpvar(nx_SE,1); % Simulated initial States x(0), this is a decision variable

u_SE = sdpvar(nu_SE, N_SE); % insulin input u(0),...,u(N_SE - 1). here they are known values.

x_hat_SE = sdpvar(nx_SE,1); % Estimated initial states x_hat(0)

y_SE = sdpvar( 1, N_SE); % measurements y(0),...,y(N_SE - 1)

% system dynamics

E = [1 0 0 0 0 0];%times xk to get the first parameter

Y_SE = [];  % simulated output H(x(t))

V_SE = [];  % V: discrepency between measurements and simulated outputs

xk_SE = x_SE;

for k = 1:N_SE

    xk_SE =  A*xk_SE + b + C*u_SE(k)  + D*w_SE(k);

    Y_SE = [Y_SE; E*xk_SE];

    V_SE = [V_SE; y_SE(k) - E*xk_SE];

end

Constraints =[Constraints,distur_1 + distur_2 + distur_3 <= 1 ];

C1 = consequtiveON([distur_1],5);

Constraints = [Constraints, C1, sum(max(diff([0;distur_1]),0)) <= 1, sum(distur_1) <= 20];

C2 = consequtiveON([distur_2],5);

Constraints = [Constraints, C2, sum(max(diff([0;distur_2]),0)) <= 1, sum(distur_2) <= 20];

C3 = consequtiveON([distur_3],5);

Constraints = [Constraints, C3, sum(max(diff([0;distur_3]),0)) <= 1, sum(distur_3) <= 20];

upperbound = sdpvar(N_SE,nw_SE);

lowerbound = sdpvar(N_SE,nw_SE);

for k = 1:N_SE

    Constraints = [Constraints, lowerbound(k) <= w_SE(k)<= upperbound(k)];

end

objective =  norm(x_SE - x_hat_SE, 1) + norm(V_SE, 1); 

% use robustmodel to pre-compile a model without solving it

ops = sdpsettings('solver','gurobi','verbose',0);

estimator = optimizer(Constraints, objective, ops, {u_SE, y_SE, x_hat_SE,upperbound,lowerbound}, {xk_SE, w_SE,objective} );

The output gave the warning: You have unbounded variables in IMPLIES leading to a lousy big-M relaxation.  What is the problem with the implies?

Thank you

[Hongkai Chen]

Then the solver became much slower.

[Johan Löfberg]

All variables in an implies operator must be explicitly bounded. You don't have bounds on distur_1_amp, w1 etc

[Johan Löfberg]

You have nothing to compare with as You've never solved a working code before

And worst-case, this model will take inf years to solve. It is a massive MILP

[Johan Löfberg]

It is crazy to try to debug etc with an horizon of 150 which leads to a huuuuge model. Start with, say, 5... Make it work for trivial cases and go from there to find the practical limit of performance

[Johan Löfberg]

this is weird as you only touch the first column

[Hongkai Chen]

Do I have to bound binvar too?

[Hongkai Chen]

upperbound, lowerbound and w_SE are all column vectors. So I guess I constrained the disturbance at all the points? Did I do it wrong? 

[Johan Löfberg]

no

[Johan Löfberg]

There are N_SE*nw_SE of them, but you only bound N_SE of them