NaN in code. Can't get around it.

[Asher Metzger]

Hi Johan,

I am grateful for you for Yalmip. I'm taking my first footstep in optimisation and Yalmip is great help.

I'm trying to model a water flow network with different constraints.

I have difficulty with the head constraint which is dependant on sdpvar decision variable. I tried the tricks you suggest in:

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.NANInModel

But haven't had success yet.

I'll be very thankful if you could offer me anther way to look at solving the NaN problem.

The head code is in line 101.

Thanks,

Asher

[Johan Löfberg]

You create h as a double, and then you assign sdpvars into it. That is precisely one of the cases discussed on the page you linked to (the loop is weird though, you loop over t but compare i with 1)

h
for t=1:2
    if i==1
        h(t+1) = ((recharge(t) - Q(1,t))/(SA)) + h_init;
    else
        h(t+1) = ((recharge(t) - Q(1,t))/(SA)) + h(t);
    end
end
h
h =

    11     0     0


h =

   NaN     0     0

BTW, why not vectorize the code?

t = 1:2
h = [h_init (recharge(t)'-Q(1,t))/(SA)]

A general comment on the model is that you introduce a nasty division with a decision variable

C0(1,t) == (1/h(t+1))*(...

You are better of with a bilinear equality instead (no risk of divide by zero when the solver works on it)

C0(1,t)*h(t+1) == (...

[Johan Löfberg]

Even weirder issue in your code is the fact that you use the variables h which simply is zero, when you define the cost function. h is not created, with actually content, until after the loop which generates the cost.

[Asher Metzger]

Thanks Johan,

I did notice that after I vectorized the h variable. I moved it up the code and it actually seems to work.

I also tried to use איק bilinear equality but I get:

Error using sdpvar/eq (line 15)

Error using constraint (line 49)

Equality constraint evaluated to trivial false (no decision variable in constraint)

any ideas?

[Johan Löfberg]

Really? I get trivial infeasibility for your original model (with h defined directly)

[Johan Löfberg]

Already these two constraints are infeasible

 optimize([C9,C10])

[Johan Löfberg]

Sorry, missed the fact that h was recursively defined. Hence, a reasonable model would be

h = h_init;
for t=1:2    
     h = [h ((recharge(t) - Q(1,t))/(SA)) + h(t)];        
end

I cannot solve the problem in that for, fmincon crashes (with divide by zero most likely)

When I use bilinear form instead, YALMIP discovers that the model is trivially infeasible (as you see also). The reason is that the expression

 C0(1,t)*h(t+1) - ((Cr*recharge(t)-Ca_init*Q(1,t))/SA + Ca_init*h(t))

evaluates to -1000, i.e, the left-hand and right-handside can impossibly be the same

[Asher Metzger]

Hi Johan,

Thanks for your answers. I realised I have to work much more organised otherwise I loose the logic of the code.

I worked on the code with my project mentor and we cleaned up the code and got it running. We still haven't found out why we get NaN's on certain circumstances. 

In the attached files we have one code which is feasible (though not yet with the wanted results) and the other getting NaN's.

In addition I get the variable Ca equal for both numbers in the vector and I think that is my problem but i don't understand why. (project_test_3... .m)

[Johan Löfberg]

Project_test_3 runs on my machine, no infeasibility, albeit many iterations
  630    2992    7.385967e+01    2.000e-08    3.405e-03    7.584e-04

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
  631    2997    7.385967e+01    1.995e-08    3.405e-03    7.583e-04

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 3000 (the default value).



test test_4, the problem is that fmincon fails to create an initial point where things can be evaluated. The problem in the model is the divisions with h(t) when creating Ca
sol=optimize(C2,[],sdpsettings('debug',1))
Error using barrier (line 26)
Nonlinear constraint function is undefined at initial point. Fmincon
cannot continue.

Error in fmincon (line 841)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] =
    barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...


You have to create a model which avoids such a nasty singularity-inducing nonlinearity at 0 (let Ca be a decision variable and multiply both sides by h(t) and introduce an equality constraint instead, etc)

[Asher Metzger]

In the latest version (attached) I avoid the singularity issues.

I still get suboptimal results. i.e. the Ca constraints don't reach max (they should).

I would like to try a different solver. what solver can replace the fmincon solver?

[Johan Löfberg]

first you have to look over the scaling of the problem. The fact that you can get very different optimal objectives by changing the magic scaling 10^6 in the objective function reveals that the model is badly formulated. After you scale it, some coefficients in the objective will be around 1, and some will be around 10^-7. This is something most solvers will struggle with (10^-7 is considered as noise by most solvers)

Are you mixing "cost of a coca-cola" and "budget of building a new wal-mart" in your objective? You are most likely using the wrong units in some part of your model. If the cost is priceofcoke*numberofcokes + costtobuildanewshop, you should probably work with the units priceofatruckofcoke*numberoftrucks + costtobuildanewshop

[Johan Löfberg]

The claim that Ca should be at its upper limit is not correct from what I can see.

If I solve the problem using fmincon, I get a solution with objective 65.7979.

If I fix Ca to Ca_max, and use a global solver (e.g., baron), it finds that the optimal objective is between 67.66 and 67.678, i.e., the lower bound from the global solver is worse than the solution in which Ca was allowed to be non-tight

Either your intuition is wrong, or the model is wrong (or the model is so horribly badly scaled that no solver can work with it robustly)

[Johan Löfberg]

and there is still a lot of weird code

On line 50 you define Cd. Then you define it again on line 64, and then you finally do something on line line 66.

Remove it all and simply use

Cd = Csea*((100-RR(t))/100);

[Johan Löfberg]

and there is a weird lack of symmetry in your code. Why isn't Ca and h handled in the same way in terms of time-indexing


Your h is simply

h = h_init
for i = 1:2
   h = [h (recharge(i)-Qa(i))/(SA)+h(i)]
end

and with a sensible indexing strategy, the Ca*h equations should be possible to write with resorting to if-cases

for t=1:2
 C5 = [C5 Ca(t+1)*h(t+1) == (Cr*recharge(t)-Ca(t)*Qa(t))/SA + Ca(t)*h(t)]; %aquifer salinity
end

[Johan Löfberg]

Additionally, since you are manually defining all those A0, Aeq, B0 matrices, you seem to hide some very simple structures in the model

The whole B0*C==0 thingy

%dilution Rule
B0=zeros(4,8);
M=[3 4; 5 6;7 8; 9 10];
for i=1:4
    B0(i,M(i,:))=[1 -1];
end
B0=[B0 zeros(4,2)];

simply says that Cp has a block structure of repeated elements. Better to define that directly and reduce the number of variables, and if you go for a global solver, reduce the number of nonlinear terms. One equivalent way is (but your particular structure can of course be done in many ways)

Cp = kron(sdpvar(4,2),[1;1])

[Asher Metzger]

Hey Johan,

Thanks a lot.

I implemented your suggestions. now the model gives an even higher result but the code is reduced significantly and I can't see anything wrong with the model.

The interesting thing is if I give the RR constraint an upper bound of 99.43 and RR(1) = 99.25 I get the anticipated result: 61.

But when there is a bigger bound range for RR [99 99.9] then I get a worse result: 82. How is this possible?

The scale is now 1:1. 

Asher

[Johan Löfberg]

nonconvexities and local minimizers

[Asher Metzger]

how then can I make sure that what I get is the absolute minimum point?

or how can I find more minimum points?

[Johan Löfberg]

Global optimization is simply ridiculously hard, so you should not expect that you can solve this problem globally, at least not without investing significant effort in various tricks to improve the model (nothing obvious as I can see) and using as much information as you possibly can to tighten bounds etc.

YALMIP built-in global solver quickly tightens the gap to say that the global solution has an objective between 46 and 66, but then it seems to stall (only ran for a minute or so though). The state-of-the-art global solver baron achives bounds of 51 and 66, but then also seems to stall.

>> sol=optimize([const,-1e4 <= recover(depends(const)) <= 1e4],sumcost,sdpsettings('solver','bmibnb','bmibnb.uppersolver','ipopt'))
* Starting YALMIP global branch & bound.
* Branch-variables : 34
* Upper solver     : ipopt
* Lower solver     : CPLEX
* LP solver        : CPLEX
 Node       Upper      Gap(%)       Lower    Open
Restoration phase is called at point that is almost feasible,
  with constraint violation 1.476986e-010. Abort.
    1 :    7.040E+01    44.08      4.461E+01   2  Improved solution  
    2 :    7.040E+01    44.08      4.461E+01   1  Infeasible  
    3 :    7.040E+01    44.08      4.461E+01   2    
    4 :    7.040E+01    44.08      4.461E+01   3    
    5 :    7.040E+01    44.08      4.461E+01   4    
    6 :    7.040E+01    44.08      4.461E+01   5    
    7 :    7.040E+01    44.08      4.461E+01   6    
    8 :    7.040E+01    44.08      4.461E+01   7    
    9 :    7.040E+01    44.08      4.461E+01   8    
   10 :    7.040E+01    44.08      4.461E+01   9    
   11 :    7.040E+01    44.08      4.461E+01  10    
   12 :    7.040E+01    44.08      4.461E+01  11    
   13 :    7.040E+01    44.08      4.461E+01  12    
   14 :    7.040E+01    44.08      4.461E+01  13    
   15 :    7.040E+01    44.08      4.461E+01  14    
   16 :    7.040E+01    44.08      4.461E+01  15    
   17 :    7.040E+01    44.08      4.461E+01  14  Infeasible  
   18 :    7.006E+01    43.63      4.461E+01  15  Improved solution  
   19 :    7.006E+01    43.63      4.461E+01  16    
   20 :    7.006E+01    43.63      4.461E+01  17    
   21 :    7.006E+01    43.63      4.461E+01  18    
   22 :    7.006E+01    43.63      4.461E+01  19    
   23 :    7.006E+01    43.63      4.461E+01  20    
   24 :    7.006E+01    43.63      4.461E+01  21    
   25 :    7.006E+01    43.63      4.461E+01  22    
   26 :    7.006E+01    43.63      4.461E+01  23    
   27 :    7.006E+01    43.63      4.461E+01  22  Infeasible  
   28 :    7.006E+01    43.63      4.461E+01  23    
   29 :    7.006E+01    43.63      4.461E+01  24    
   30 :    7.006E+01    43.63      4.461E+01  25    
   31 :    7.006E+01    43.63      4.461E+01  26    
   32 :    7.006E+01    43.63      4.461E+01  25  Infeasible  
   33 :    7.006E+01    43.63      4.461E+01  26    
   34 :    7.006E+01    43.63      4.461E+01  27    
   35 :    7.006E+01    43.63      4.461E+01  28    
   36 :    7.006E+01    43.63      4.461E+01  29    
   37 :    7.006E+01    43.63      4.461E+01  30    
   38 :    7.006E+01    43.63      4.461E+01  31    
   39 :    7.006E+01    43.63      4.461E+01  32    
   40 :    7.006E+01    43.63      4.461E+01  33    
   41 :    7.006E+01    43.63      4.461E+01  34    
   42 :    7.006E+01    43.63      4.461E+01  35    
   43 :    7.006E+01    43.63      4.461E+01  36    
   44 :    7.006E+01    43.63      4.461E+01  37    
   45 :    7.006E+01    43.63      4.461E+01  38    
   46 :    7.006E+01    43.63      4.461E+01  39    
   47 :    7.006E+01    43.63      4.461E+01  40    
   48 :    7.006E+01    43.63      4.461E+01  41    
   49 :    7.006E+01    43.63      4.461E+01  42    
   50 :    7.006E+01    43.63      4.461E+01  43    
   51 :    7.001E+01    43.55      4.461E+01  44  Improved solution  
   52 :    7.001E+01    43.55      4.461E+01  45    
   53 :    7.001E+01    43.55      4.461E+01  46    
   54 :    7.001E+01    43.55      4.461E+01  47    
   55 :    7.001E+01    43.55      4.461E+01  48    
   56 :    6.762E+01    40.28      4.461E+01  49  Improved solution  
   57 :    6.762E+01    40.28      4.461E+01  50    
   58 :    6.762E+01    40.28      4.461E+01  51    
   59 :    6.762E+01    40.28      4.461E+01  52    
   60 :    6.762E+01    40.28      4.461E+01  53    
   61 :    6.762E+01    40.28      4.461E+01  54    
   62 :    6.762E+01    40.28      4.461E+01  55    
   63 :    6.762E+01    40.25      4.463E+01  56    
   64 :    6.762E+01    40.25      4.463E+01  57    
   65 :    6.762E+01    40.25      4.463E+01  58    
   66 :    6.762E+01    40.25      4.463E+01  59    
   67 :    6.762E+01    40.25      4.463E+01  60    
   68 :    6.762E+01    40.25      4.463E+01  61    
   69 :    6.762E+01    40.25      4.463E+01  62    
   70 :    6.762E+01    40.25      4.463E+01  63    
   71 :    6.762E+01    40.25      4.463E+01  64    
   72 :    6.762E+01    40.25      4.463E+01  65    
   73 :    6.762E+01    40.25      4.463E+01  66    
   74 :    6.762E+01    40.25      4.463E+01  67    
   75 :    6.762E+01    40.25      4.463E+01  68    
   76 :    6.675E+01    39.03      4.463E+01  69  Improved solution  
   77 :    6.675E+01    39.03      4.463E+01  70    
   78 :    6.675E+01    39.03      4.463E+01  71    
   79 :    6.675E+01    39.03      4.463E+01  72    
   80 :    6.675E+01    39.03      4.463E+01  73    
   81 :    6.675E+01    39.03      4.463E+01  74    
   82 :    6.675E+01    39.03      4.463E+01  75    
   83 :    6.675E+01    39.03      4.463E+01  76    
   84 :    6.675E+01    39.03      4.463E+01  77    
   85 :    6.675E+01    39.03      4.463E+01  78    
   86 :    6.675E+01    39.03      4.463E+01  79    
   87 :    6.675E+01    39.02      4.463E+01  80    
   88 :    6.675E+01    39.02      4.463E+01  81    
   89 :    6.675E+01    39.02      4.463E+01  82    
   90 :    6.675E+01    39.02      4.463E+01  83    
   91 :    6.675E+01    39.02      4.463E+01  84    
   92 :    6.675E+01    39.02      4.463E+01  85    
   93 :    6.675E+01    39.02      4.463E+01  86    
   94 :    6.675E+01    39.02      4.463E+01  87    
   95 :    6.675E+01    39.02      4.463E+01  88    
   96 :    6.675E+01    39.02      4.463E+01  89    
   97 :    6.675E+01    39.02      4.463E+01  90    
   98 :    6.675E+01    39.02      4.463E+01  91    
   99 :    6.675E+01    39.02      4.463E+01  92    
  100 :    6.675E+01    39.02      4.463E+01  93    
* Finished.  Cost: 66.7534 Gap: 39.0235
* Timing: 23% spent in upper solver (96 problems solved)
*         2% spent in lower solver (164 problems solved)
*         67% spent in LP-based domain reduction (6242 problems solved)