Cost function minimization wrt one decision variable

[Dario Zurlo]

Hello,

I am trying to solve the optimization problem in the screenshot. My problem is

that I want to minimize the cost function wrt p_i but I think that if I do not specify

anything the minimization is wrt all the decision variables, that in my case are v and p.

I tried to set v as parameter when I call the optimizer but seems it does not work.

It is possible to specify which decision variable I want to use for the minimization?

Or I have to define v as another type?

Regards

Zurlo

[Johan Löfberg]

What does it even mean to not minimize over a subset of the variables? If the solver cannot give them values, do you mean they are constants (why even talk about them then), or forced to take on a fixed value (then you model is missing some equalities), or uncertain and you want to solve a minmax problem (does not look like that)

If you have a truck and the decision variables are size of engine and length of truck (which to an extent is a function of the engine size as a larger engine requires a larger space for the engine, but you can also design the loading area) , and your objective is to minimize prize. Now you say that you only want to minimize over engine size, and you are not allowed to actively minimize over he other variables? The cheapest solution is the truck with the smallest engine and as short as possible, but you want to find a solution where we only optimize over engine size and the length of the car just magically isn't selected (what value should it have then? if we cannot select it in the optimization the cost cannot even be defined . Is it fixed a-priori but then it is just a constant and shouldn't be called a decision variable )

To me, your problem looks like a nonlinear nonconvex MPC problem with control input pk and state-sequence vk. You should remove the possible singularity by multiplying the state-update with vk to arrive at a bilinear equation instead, and write the nonconvex nonsmooth non-zero condition as vk^2>=e^2, leading to a smooth nonconvex quadratically constrained problem instead.

[Dario Zurlo]

Thanks for the quick response.

Yes exactly my problem is a nonlinear optimization problem, in which the first constraint is the discrete dynamical model, discretized wrt the space Deltax and the second constraint needs for avoid singularity.

The input is pk and the states are xk and vk, I have to implement an MPC but I'm trying to solve just in one iteration. Anyway I think about what you told me and I put this optimization problem in a MPC scheme.

Thanks

[Johan Löfberg]

if you have both the states and inputs as decision variables (i.e. both sdpvars and defining the dynamics using equality constraints) then they are, by definition, both decision variables. If you have only defined the input as decision variable and defined states using assignments, well then only the input is a decision variable as everything else is expressed in the input

% Only one decision variable x, otimization over x

sdpvar x

y = x + 1;

optimize(x<= 5,x^2)

% Two decision variables, optimization over x and y

sdpvar x y

optimize([y == x + 1;x<= 5],x^2)

[Johan Löfberg]

was supposed to be

% Only one decision variable x, otimization over x

sdpvar x

y = x + 1;

optimize(y<= 5,x^2)

% Two decision variables, optimization over x and y

sdpvar x y

optimize([y == x + 1;y<= 5],x^2)

[Dario Zurlo]

This is my Matlab code

%Parameters
beta = 50;
g = 9.81;
z = 0.8;
w = sqrt(g/z);
Deltax = 0.001;
epsilon = 0.005;
P = 202;


%decision variables
%x=sdpvar(2*ones(1,P+1),ones(1,P+1));
u=sdpvar(1*ones(1,P+1),ones(1,P+1));


%To define constraints and the cost function I iterate over all the control horizon
Constraints = [];
for i = 1:P
        Constraints=[Constraints;
        x{i+1}(2)*x{i}(2) == x{i}(2)*x{i}(2) + w.^2*(x{i}(1)-u{i}(1)) * Deltax;x{i}(2).^2>=epsilon.^2;];
end
Cost=0;
for i=1:P
    Cost = Cost + (x{i}(2)-vref(i)).^2+beta*(u{i}-pref(i)).^2;
end

Problem = optimize(Constraints, Cost);

I tried both x as decision variable, but I think that is wrong because it is not a decision variable but only a state, and as it is in the code above. In both cases I obtain the same result

[Johan Löfberg]

x is a decision variable. That does not mean it can be selected arbitrarily as it has to satisfy the state-dynamics

Just like this example, it might be that x is the actual decision we make (the number of items we buy, we have one item in storage right now so next day we will have x+1 items), and y is a direct consequence of that choice. That does not prevent us from defining y as a decision variable, the  solver has no freedom but forcing it to be x+1. 

sdpvar x y

optimize([y == x + 1;y<= 5],x^2)

Looks like you have manually discretized a nonlinear optimal control problem. You are probably much better of using a dedicated optimal control software package, which typically will do much more clever discretization/parameterization/structure exploitation etc.

[Johan Löfberg]

you say "I tried both x as decision variable, but I think that is wrong because it is not a decision variable but only a state, and as it is in the code above". That is impossible as the code as shown does not run as x isn't defined. If it ran, then just just happened to have the variable still defined in the work-space, meaning you ran the same thing again

fredag 5 mars 2021 kl. 11:41:00 UTC+1 skrev zurlo...@gmail.com:

[Dario Zurlo]

yes x is in the work space for this reason "works" I forgot to clear the workspace. Then if I don't want to define x as decision variable I need to impose x as the dynamic equation and also impose in the constrint, isn't it?, but as you told me even if x is a decision variable it cannot be arbitrary, of course because it is constrained by the dynamic.

I have discretized the system manually because I want to perform a spatial quantization and not a time quantization as usual as it is possible to see from the figure, for this reason the time step is not constant and it is inversely proportional to the state vk.

thanks for your patience

[Johan Löfberg]

First, your problem is not well-posed as you are missing an initial condition on x, i.e. when optimizing over both x and u either x{1}==something or simply setting it to a value before definng the state-dynaics constraints, x{1}=something

Also, the model you are setting up make no sense, as you are only defining the dynamics of the first state in the vector x{i}. For that reason, I am going to assume it really is supposed to be a scalar state

If you would want to setup the whole model without having x as a decision variable (which you don't, a explicit representation in control input only is typically a bad numerical model, already in the linear case, much better to create a sparse model in state and control input as they typically have better numerical properties, and in the nonlinear case they simply don't work, as can be seen below) you cannot work with the bilinearized model, as it is based on an implicit representation through products of current and next state. Hence you have to keep the original nasty division

Iterate through the dynamics with a model which expresses everything in the control space only. Note that this expression simply becomes exceedingly complex, and it leads to a model which simply cannot be used. That's why I am printing the expression so you understands why it stalls (the code below will not terminate, your computer will run out of memory). For instance x{21} is an expression which contains close to 13 million characters (which indicates how many monomials etc there are in that term). That's why you use an implicit representation with both x and u as decision variables

x{1} = 1; % Or whatever you have

for i = 1:P

x{i+1} = x{i} + w.^2*(x{i}(1)-u{i}(1))*Deltax/x{i};

expr = sdisplay(x{i+1});

    expr{1}

    length(expr{1})

end

[Johan Löfberg]

and what you don't even see is that even if you think you are creating a model which only lives in the u-space, under the hood the division can introduce new internal variables and constraints that would be added to the dynamic model. This happens when the denominator is anything but a simply monomial

if you for instance create something like this (which is the kind of stuff you will have when the iteration starts moving forward)

x(2) = u(2)/(1 + u(1))

YALMIP will internally flatten this to x(2) = u(2)*z and then internally remember the constraint z*(1 + u(1)) == 1 which will be added to the model. I.e. even if you think you have only have the u(1) and u(2) here, it would be u(1), u(2) and z 

[Dario Zurlo]

At the end I have defined all the variables as decision variables, because there is no reason not to as you made my understand. I have also defined the constraint for the first state that is x{i+1}(1) == x{i}(1)+Deltax; and of course I have defined the initial state and the initial input that are x{1}(1) = 0; x{1}(2) = epsilon; u{1} = 0;. The result now seems better even if the vk profile is not as expected but this is another problem, may be I have to tune the weight beta in the cost function or the spatial constant, instead the xk profile seems good.

I have addressed the same problem quantized  wrt the time and quantized wrt the space the latter has caused me a lot of confusion, thanks for the patience.

[Johan Löfberg]

You have to have in mind that the solver could have failed completely and returned something far from optimal

You should start with very short horizon and compare the solution with the globally optimal solution computed using gurobi or bmibnb etc, to see if your local solutions are in the right ballpark. If they are way off already for a trivial prediction length, you cannot have much confidence when going up

[Dario Zurlo]

it was just a stupid mistake now all the result are as expected, thanks for the advice you gave me. This the code:

beta = 1.50;

g = 9.81;
z = 0.8;
w = sqrt(g/z);
Deltax = 0.001;
epsilon = 0.005;

P = 1500;

x=sdpvar(2*ones(1,P+1),ones(1,P+1));
u=sdpvar(ones(1,P+1),ones(1,P+1));

x{1}(1) = 0;
x{1}(2) = epsilon;
u{1} = 0;

Constraints = [];
for i = 1:P

        Constraints=[Constraints;

            x{i+1}(1) == x{i}(1)+Deltax;

            x{i+1}(2) == x{i}(2) + w.^2*(x{i}(1)-u{i}) * Deltax/x{i}(2);

            x{i}(2).^2>=epsilon.^2;];


end

    Constraints=[Constraints;x{P+1}(1) == pref(1,P+1)];
Cost=0;
for i=1:P
    Cost = Cost + (x{i}(2)-vref(1,i)).^2+beta*(u{i}-pref(1,i)).^2;
end

options = sdpsettings('verbose',1);

Problem = optimize(Constraints, Cost,options);

The only problem, unavoidable, is the execution time because the horizon is very big P = 1500

[Johan Löfberg]

why are you complicating matters with the first state? it's simply (i-1)*Deltax

[Johan Löfberg]

and as expected, testing with a short horizon P=10 and just some random vref pref which you haven't supplied shows that a local solver (fmincon, ippt) easily fails to find the global solution and delivers a sub-optimal one

Also, if you really want to use long horizons, you should really vectorize the code. It is messy to do when you are using cells, and there is no real reason to use cells here, so if you simply use standard vectors instead, you can write the state dynamics as (dumping the first state which is predefined, and bilinearizing it to get rid of the singularity inducing division)

x(2:end).*x(1:end-1) == x(1:end-1).^2+w.^2*((0:P-1)*Deltax-u(1:P))*Deltax

similarly the objective is a simply sum(().^2)

I would also try to dump the x^2>=eps^2 constraints. It appears to be satisfied anyway, at least solve first and then add if necessary

lördag 6 mars 2021 kl. 18:13:34 UTC+1 skrev zurlo...@gmail.com:

[Dario Zurlo]

 the position x at step i is the discretization step  Deltax times the step, then xi = Deltax*i no (i-1) that is simpler of I had written in the code, now I change it.

Do you mean that the solution that I have found may be no the optimal one? Even if it is no the optimal solution I obtain a very good solution that is satisfactory for me.

Which solver I should use?

Now I try to modify also the code in the vectorized form that looks neater.

[Johan Löfberg]

Since you are using a local solver I guess (fmincon, ipopt) you have no guarantees what so ever. I would at least try and compare with a global solver (baron, bmibnb, gurobi) to see how well the local solver performs for shorter horizons