Counting an occurrence inside the optimization

[Louis Hallerdt]

Hello,

I have an sdpvar "socdv" in my code and would like to know if it is possible to count inside the optimization, the number of time its derivative socdv(:,t+1) - socdv(:,t) becomes negative from a local maximum value of socdv at time t.
Is it possible to have a binary variable that becomes true every time the conditions mentionned above are satisfied.
Outside the optimization it is easy since I can use "if/else" command but with sdpvar that do not have any value before the solving, this method is not possible.

thank you in advance.

[Johan Löfberg]

It is unclear what you want to do. Show how you would define this value if this variable was a simple double

[Louis Hallerdt]

"SOC" is the double equivalent of "socdv" which are of size (Ns,timestep+1)

In double, what I would like to do is: Counting the number of discharge cycles with the variable "Ncycle"

I attached an example of the code.

[Louis Hallerdt]

The link didn't work, here is the code

[Louis Hallerdt]

After having found this functionality on the command "implies" (http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.Implies), I tried to model the cycle counter like this as a constraint:

cycleornot = binvar(Ns,timestep);
for b=1:Ns
    for t=1:timestep
        CycleCONDITION = [socdv(b,t+1) - socdv(b,t) <= -10^(-4)];
        F = implies(CycleCONDITION, cycleornot(b,t));
        Constraints = [Constraints, F];
    end
end

I expected this code to give me, after optimization, true value for battery "b" at time "t" when a cycle begins.
As the binary variable "cycleornot" does not affect the cost in any way (and nor the other constraints than "F"), I thought I would get the same optimal cost and "socdv" values as the initial code without the few lines above.
Now my results are wrong and the optimizer fails.

I used these settings to solve my problem:

options = sdpsettings('solver','gurobi','gurobi.QCPDual',1);
sol = solvesdp(Constraints,Objective,options);

Where have I done something wrong?
Thank you in advance for helping me.

[Johan Löfberg]

You are correct that you must use logic programming with implies/iff etc to model this.

The fact that it oes infeasible when you add these activvation constraints is weird. I can see two reasons

1. Ns=timestep so you should use 'full' to ensure that cycleornot is fully parameterized. Still don't see how that could possibly affect feasibility, as all ones is feasile

2. Most likely: You problem is badly sclaed and the optimal socdv has elements larger than 10^4, which is the default bound added by YALMIP when it performs big-M modlling, and you haven't explicitly added any variable bounds on the variables involved in the implies operator

If you post reproducible code, it might be possible to see the issue

[Louis Hallerdt]

I cleared Matlab workspace then the objective function got the same result as before the usage of the binary variable.
Here is a fully reproducible code which only has an issue on counting the number of battery cycle using the binary variable "cycleornot"

Thank you very much in advance for helping me

[Johan Löfberg]

I cannot run the code (cannot perform && to concatenate constraints on line 289) so I am curious in which YALMIP version / MATLAB version this actually runs?

Anyway, edited it by replacing && with , and it runs. However, it is feasible, which you seem to say in the file also, so I don't see what the problem is?

BTW, this

Constraints = [Constraints, cycleornot(:,1)==repmat(0,Ns,1)];

is more conveniently written as

Constraints = [Constraints, cycleornot(:,1)==0];

and solvesdp is obsolete, it is called optimize now, and double is called value.

[Louis Hallerdt]

Good evening Johan,

thank you for your fast reply.
I installed YALMIP in February, I don't think my version is too outdated. I use MATLAB R2014b.

the problem is the "cycleornot" is set to 1 way too often.
In the current state of the program, it should be 1 for only two timesteps for each battery.

A cycle is modelled to start whenever the derivative of the SOC becomes negative from a local maximum
When you look at the top-left plots on Figure 11 and 22, the SOC satisfies these conditions at time 10 and 17 for both batteries so "Ncycles" should be [2;2] instead of [23;23] ...

i must have given uncomplete or wrong conditions in implies at line 298.

[Johan Löfberg]

cycleornot all ones is feasible regardless of the soc variables. You seem to confuse implication with if-and-only-if (iff).

However, I would reverse the logic, since these things are terribly sensitive to tolerances in solvers. implies(x>=0,d) might return x = 0.000000001 and d=0, as that is feasible up to solver tolerances. Typically, it is better to try to use implies(binary,condition), and then explicitly model disjunct statements using explicit disjunctions such as implies(d1, x>=0)+implies(d2,x<=0)+(d1+d2==1), as you then have a clear definition of what the solver actually means with the solution. Finally, you have to have a clear idea of wether you will use the solution which was computed to a-posterio define the values of the logic values, or if you will use the binary values defined by the solver (which might suffer from tolerance issues)

[Louis Hallerdt]

I checked the iff command and that is exactly what I wanted to tell the solver.
So this line means cycleornot(b,t) will be feasible (true) when the constraint vector [CycleCONDITION1, CycleCONDITION2] is satisfied, isn't it?

F = iff([CycleCONDITION1, CycleCONDITION2], cycleornot(b,t));

But now cycleornot(b,t) only became true at b=2, t=10 whereas I want cycleornot(1,10),cycleornot(1,17) and cycleornot(2,17) to be also 1.
I intend to use the results from cycleornot for further optimization of the objective which will have a cost depending on this number of cycles this is why I am trying to count the cycles inside the optimization.

The binary variable is, at the moment, only used to get the information of number of cycles and my constraints are not piecewise defined so the implies command may not be usefull in this way.

[Johan Löfberg]

this is precisely he issue I am talking about. It is extremely hard to activate a binary variable based on a condition of a continuous variable, in particular if the variable can be big (thus inducing a large big-M constant). Change the magic number from 10^-4 to 10^-2 and you'll see that it all of a sudden behaves as expected

also note that gurobi actually warns about significant constraint violations

[Louis Hallerdt]

It worked! However the computing time has increased by a factor of 10 (about 18seconds now vs few seconds before)
Is it due to the complexity of operations the binary variables and iff involve?
When you say "extremely hard to activate a binary variable based on a condition of a continuous variable' for the solver, is there other consequences else than the longer computation time?

Below is the message I get from the solving

Optimize a model with 14982 rows, 11852 columns and 33568 nonzeros
Model has 888 quadratic constraints
Coefficient statistics:
  Matrix range    [2e-08, 3e+02]
  Objective range [6e-02, 1e+05]
  Bounds range    [1e+00, 1e+00]
  RHS range       [2e-04, 1e+02]
Presolve removed 6648 rows and 2022 columns
Presolve time: 0.29s
Presolved: 11838 rows, 13334 columns, 32662 nonzeros
Variable types: 13150 continuous, 184 integer (184 binary)

Root relaxation: objective -5.174086e+03, 8422 iterations, 0.92 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

H    0     0                    81560.762881          -      -     -    3s
     0     0 3765.37827    0    6 81560.7629 3765.37827  95.4%     -    4s
     0     0 34666.4068    0   27 81560.7629 34666.4068  57.5%     -   10s
     0     0 35308.7177    0   10 81560.7629 35308.7177  56.7%     -   10s
     0     0 35308.7177    0    6 81560.7629 35308.7177  56.7%     -   11s
     0     2 35308.7177    0    6 81560.7629 35308.7177  56.7%     -   11s
*   18    10               6    35754.230725 35714.2777  0.11%   174   14s
    22     4     cutoff    6      35754.2307 35714.2777  0.11%   179   15s

Cutting planes:
  Learned: 4
  Clique: 2

Explored 29 nodes (19221 simplex iterations) in 18.36 seconds
Thread count was 2 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Warning: max constraint violation (2.0055e-05) exceeds tolerance
Best objective 3.575423072471e+04, best bound 3.575423072471e+04, gap 0.0%

[Johan Löfberg]

when i say extremely hard I mean from a numerical stability point of view. Of course iff adds huge numbers of binaries, so that makes the problem hard from a complexity point of view too

[Louis Hallerdt]

I understand.
Tusen tack for all your help and advices!