Re: optimizer in parametric problem using YALMIP

[Johan Löfberg]

Code does not run (missing many parameter values)

feels like you are sort of re-inventing the wheel though. Why not simply solve the nonlinear program (either using a local or global solver)

On Monday 2 September 2024 at 12:27:20 UTC+2 georgia...@gmail.com wrote:

Hello!

I am quite new to the optimization world, as well as Matlab-YALMIP and I am trying to model a parametric maximization problem using YALMIP, which is solved iteratively within an iterative convexification method. However, Matlab crushes even for a very small instance of the problem (N = 3, timesteps = 10,  N_scen = 3).  Also, for the cases I manage to get an output, each iteration is slow. My parameters are C_max, B, P_h, rho_scen, h_scen, s_j, C and my optimization variables r, l, p, SoE, P_a. I use piecewise linearization to construct logarithms of variables and of expressions containing parameters and variables both in my objective function and constraints. How could I speed up the performance of optimizer?  Do you think that there are alternative ways to express the logarithms so to get a solution faster?

Thank you very much, in advance!!

[Johan Löfberg]

If you profile the code you will see that it stalls due to poor vectorization

Massive time is spent in this trivial elementwise loop

for m = 1:M

approx_constraints = [approx_constraints, segment_values(m) == c(m) * x_var + h(m)];

end

which simply should be written using segment_values == c* x_var + h

and similiar loops below

On Monday 2 September 2024 at 14:05:43 UTC+2 georgia...@gmail.com wrote:

I want to compare the iterative method which solves a convex problem (with gurobi) at each iteration with the performance of a solver that solves directly the non-convex problem. I am adding the code for the method (which sets the parameter values for the toy example) that solves iteratively the problem formulated in prob.m. I would like to solve much larger instances of the problem, that is why I am trying Yalmip and optimizer to avoid defining the problem repeatedly and speed up the solution at each iteration and when the values of parameters change.

[Georgia S11]

Thank you for the advice, it sped up the problem definition (line where i just call prob), however each iteration (contb in iter.m takes too much time) is really slow and matlab crushes again. I should expect that each iteration runs really fast because the problem has been defined (where prob is called) and the only difference is that parameter C changes at each iteration. Have i misunderstood the way optimizer works? The corresponding code in python (at each iteration the problem is defined again, parameters are constants (not vars) and the problem is solved from scratch, gurobi also does the piecewise linearization with built in functions) works much faster. 

[Johan Löfberg]

You seem to create a very complex parameterization, which thus will be very expensive to construct and reason with

Also, I don't see how this can become a linear model.

You create this

expr = 1 + (h_scen(i, j, t, s) * p(i, j, t, s) * reciprocal(i, j, t, s));

where the only parameter are

{C_max; B; P_h; rho_scen; h_scen; s_j; C}

hence expr is bilinear in decision variables p and reciprocal

then you construct constraints with x_var replacing expr as

approx_constraints = [approx_constraints, segment_values == c(:) * x_var + h(:)];

so bilinear equality in p and reciprocal for fixed c

[Johan Löfberg]

...and the develop version of yalmip supports the pwa approximation stuff in gurobi 11, with log being one of the supported cases. it defaults to a true global approach but you can change the options in gurobi to used fixed intervals

On Monday 2 September 2024 at 15:40:42 UTC+2 georgia...@gmail.com wrote:

[Georgia S11]

I used this reciprocal as an auxiliary variable, because I thought optimizer treats parameters (defined with sdpvar) as variables and I wanted to avoid writing the expressions like 1+(h(i,j,t,s)*p(i,j,t,s)/(s_j + C(i,j,t,s)) (h, s_j, C are parameters and only p is the variable), which is what I want to represent . If I write it in this way I get an error "GUROBI does not support general polynomial expressions". Also in your last message, are there any built in function to model this pwa approximation?

Thanks again for the help!

[Johan Löfberg]

if you have things like x*a + b/a where x and y are decision variables, and a is a parameter, you should see it as x*a+y*b where a and b are parameters (and b happens to be a parameter which takes values tied to a)

[Johan Löfberg]

meant x*a + y/a which should be seen as x*a + y*b

[Georgia S11]

I made the change you proposed by combining the parameters of expr to one (m_scen = h_scen/s_j+C) so expr becomes 1+m_scen*p. However, each iteration is still really slow (took 40 minutes for the first iteration for N = 3, timesteps = 10, N_scen = 5, which is a really small instance). If I try to use directly log2(1+m_scen*p) i get the error "Parameters are currently only allowed to enter function such as exp, sin etc as exp(a), sin(b) etc.". I don't know which commands do the pwa directly in matlab for the logarithm (in python there is a function model.addGenConstrLogA), but apart from the variable I also have a parameter inside the logarithm, which yalmip does not seem to support and that is why I tried to express it with this custom piecewise_log2 function.

[Johan Löfberg]

I cannot debug without the code

If you simply want gurobi to use pwa approximation of log, which it supports in v 11, you can do that directly in the develop branch

Default is adaptive (i.e. vanilla global solver similiar to YALMIPs bmibnb, and solvers such as baron). If you want to use fixed intervals you have to switch to that
optimize([10 >= x >= 1], log(x),sdpsettings('solver','gurobi','gurobi.funcnonlinear',0))

[Johan Löfberg]

I've checked in two silly bottlle-necks which speeds up things for massive parameterizations like this, so use develop branch

You still have inefficient loops. for instance, sum_term2 should be done as

e = p(:,:,:,1);

sum_term2 = V*sum(e(:).^2);

With that, gurobi is now the major dominant bottle-neck (you should set verbose to 2 in creation of optimizer object, so you see when gurobi is working)

[Georgia S11]

I removed all the loops by using colons and masks to represent "if" conditions and I just used log2 directly to express logarithms (without using the funcnonlinear parameter), as I think you suggested. However,  the solution still gets a really long time to be found (less time but still impractical). In python, where I just set the parameters and each iteration re-defines and re-solves the problem (without exploiting the similarity of the problems solved in each iteration, which Yalmip and optimizer does) takes just a few seconds for an instance like this (N = 3, timesteps = 3, N_scen = 5), using inefficient loops like the ones I had in the previous messages. I wonder why this happens, because the reason I thought of using optimizer is to speed up my code even more, by formulating the problem as a parametric one and create larger instances. In the iter.m where I call prob, the problem is defined really quickly, while each iteration is performed really slow. I expected the opposite. Do you have any idea, why this happens and why the same problem in python is much faster? Do I use optimizer wrong?

[Johan Löfberg]

If you are using the gurobi framework directly, and it takes gurobi much longer to solve the problem, than when you setup the model in python and call gurobi, well then you are setting up a different problem, or you simply facing the problem that a very very small change to a MILP can make massive differences to the solution time. 

If you keep the funcnonlinear setting to YALMIPs default value (1), you are forcing gurobi to solve the problem globally, which is very dfifferent from gurobis default value (0) where it simply sets up the PWA approximation once and then solves that problem.

[Johan Löfberg]

and it goes without saying that you should of course look at the display by gurobi to see that you actually have the expected number of variables, constraints, and structure

[Georgia S11]

Thank you very much for the insights! I was mistaken in my previous message though, I called a version of a problem that used the manual piecewise approximation with the breakpoints, slopes and intervals and that's why i got a result. Using log does not work directly (for example log2(1+m_scen.*p) where p are decision variables and m_scen parameters, both set as sdpvars, I get the error "Incorrect number or types of inputs or outputs for function log".). I have an academic license of gurobi 11 (granted on February) and cloned the develop branch of yalmip one day ago. Did you mean something else that "the develop version of yalmip supports the pwa approximation stuff in gurobi 11, with log being one of the supported cases"? I thought that if I downloaded again yalmip and just directly use log2, this would work. Should I write it differently or include something more on the code?

[Georgia S11]

Hello Prof. Johan,

Thank you again for your invaluable help and the insights you’ve provided regarding my problem, YALMIP and Gurobi. I really appreciate your support.

I am really sorry for bothering you again, but I would like to ask one more time whether I can use logarithms directly on sdpvars for a problem that is defined to be solved with optimizer or I need to continue with the custom way I was piece-wise approximating logarithms, defining slopes and intercepts. 

I realised that the error I was explaining in my previous message was due to the vectorized input I used for the logarithm, so I kept the vectorized expressions only in other constraints and used for-loops for the logarithms like this one (I don't include again code for debugging, just the correspondent lines):

for i = 1:N

    for j = 1:N

        for t = 1:timesteps

            for s = 1:N_scen

                if adj_matrix(i, j) == 1

                    sum(r(i, j, :, t, s)) <= B * log2(1 + (m_scen(i, j, t, s) * p(i, j, t, s)));

                end    

            end

        end

    end

end

In the above code r, m_scen, p are defined as sdpvars (with r, p being decision variables and m_scen parameters of the problem).

From your previous messaged I understood that if I used the last versions of YALMIP and Gurobi and the parameter 'gurobi.funcnonlinear' equal to 0,  this would work as it is, but I still get the error:

“Error using optimizer (line 313)

Parameters are currently only allowed to enter function such as exp, sin etc as exp(a), sin(b) etc.Parameters are currently only allowed to enter function such as exp, sin etc as exp(a), sin(b) etc.”

Is there a point in your explanation which I have not understood? Do I need to include something more so that gurobi does directly pwa of logarithms of expressions containing sdpvars? I just want to know if this is something feasible in Matlab, Yalmip and Gurobi and if yes, how should one write such an expression so that it works.

Thank you once again for your time and guidance and really sorry for the long thread!

[Johan Löfberg]

This is not allowed

sdpvar x a
Model = log(1 + x*a) >= 1
P = optimizer(Model,x,sdpsettings('solver','gurobi'),a,x)

since there is a limitation disallowing parameters inside nonlinear functions to simplify things. Hence, you go around it by introducing a new variable

sdpvar x a y
Model = [y == x*a, log(1 + y) >= 1]
P = optimizer(Model,x,sdpsettings('solver','gurobi'),a,x)

[Georgia S11]

Thanks for the suggestion! 

I added:

y = sdpvar(N, N, T/Delta_t, N_scen, 'full');

%Rest of the code

Constraints = [Constraints, y == m_scen .* p];

for i = 1:N

for j = 1:N

for t = 1:timesteps

for s = 1:N_scen

if adj_matrix(i, j) == 1

Constraints = [Constraints, sum(r(i, j, :, t, s)) <= B * log2(1 + y(i, j, t, s))];

end

end

end

end

end

and i got again the error:

"Error using optimizer (line 313)

Parameters are currently only allowed to enter function such as exp, sin etc as exp(a), sin(b) etc."

same for a part of my objective function (where I did the same):

x = sdpvar(N, N, T/Delta_t, 'full');

Constraints = [Constraints, x == l];

logarithm = 0;

for i = 1:N

for d = 1:N

for t = 1:timesteps

logarithm = logarithm + log2(x(i, d, t));

end

end

end

Objective = logarithm + ...

(Applying log2(l(i,d,t), i.e., directly in the decision variable without using the auxiliary x, also outputs the same error. However l(:,:,:) are not parameters, but variables.)

The optimizer function is defined like this contb = optimizer(Constraints, -Objective, options, {P_h; rho_scen; m_scen}, {r; l; p; SoE; P_a; Objective});

[Johan Löfberg]

Are you really using the develop branch?

[Johan Löfberg]

Now I see the problem, appears to be something strange in the develop branch causing it to refuse

[Johan Löfberg]

See now. The support for nonlinear operators in gurobi is very recent, and sort of off-the-standard track of how solvers are setup and work with nonlinearities. At the moment, optimizer simply get confused by gurobi being used for a nonlinear operator

[Georgia S11]

So, if I understand correctly, currently the only way to define a problem including logarithms to be solved by optimizer and gurobi is to manually perform piecewise approximation of logarithms, defining slopes and intercepts, right?

[Johan Löfberg]

edit optimizer.m line 313 to "if 0" and then use a simåle variable log2(y) via an equality, and it should run

[Georgia S11]

Thanks for all the feedback you have provided, indeed MATLAB runs if I modify optimizer, but it has different results than the one I get in Python and it is also slower. I suspect that this has to do again with the logarithm. In MATLAB I use directly log2() as we said, and in Python the addGenConstrLog() function with FuncPieceError = 0.001. In MATLAB if I set FunNonlinear to 0 my method does not converge (for an instance that in Python converges in 5 seconds) but each iteration is really quick, and if I set it to 1 each iteration in slower but I get a result (more delayed than in Python). Apart from that the way I am defining the logarithm and the fact that the problem in Matlab is parametric with parameters defined as sdpvars and solved with optimizer,  it is identical in both languages. Do you think that editing optimizer and allowing to do something it is not supposed to, messed up with the solver behavior? Is there something else that I haven't consider? Should I use Python at last? I just thought that optimizer would solve my problem quicker by defining it as a parametric problem, and would allow me to solve my problem for large instances.

[Johan Löfberg]

You are comparing apples to oranges. If you set funcnonlinear to 0, you solve 1 MILP approximation only, and then it is up to you how to set up approximations. If you set funcnonlinear to 1, you are solving a global optimization problem using convex envelopes, which is much much harder than simply setting up a single approximation.

Using optimizer, which requires you to know exactly what goes on with possibility to verify that the instantiated problem really is what you want when the parameter is fixed, together with experimental stuff like the pwa stuff in gurobi, sounds like way too many layers of modelling if you really want to setup something very specific for comparison

[Georgia S11]

But if I use directly the log2() function for logarithms by editing the optimizer as we said, how can I define the piecewise approximation? Isn't the function calculating directly the logarithm? In Python I can use this addGenConstrLog() function which does a piecewise approximation with an error that can be set by the parameter FuncPieceError. How can I set up the approximation in Matlab when using the log2() function and funcnonlinear to 0?

Thanks for the help!

[Johan Löfberg]

If you have funcnonlinear set to 1 (which yalmip defaults to when creating options for gurobi) then gurobi will implement a global spatial branch-and-bound strategy like this https://yalmip.github.io/tutorial/envelopesinbmibnb. If you set it to 0 (which is normal default in gurobi) then gurobi will simply create a pwa approximation of the log, with precision controlled via the various options available for this

[Georgia S11]

Thank you for the clarification! If I remove the logarithm my program in the two languages has the same output which means that it is defined in the same way, however when I add the logarithm, with the same parameters, the result is different. Is it possible that gurobi handles logarithms differently in the two languages or yalmip calls gurobi differently? I can't think of anything else.

[Johan Löfberg]

you will have to put breaks in code so you can see the actual call being made to gurobi. There are so many things that can differ, as YALMIP introduces various auxilliary variables to simplify model generation. Add to that the complication of using an optimizer construct and a very recent non-standard add-on to gurobi