Unexepected implies warning: cannot resolve warning and unexpected solution

[mrkento11]

Dear Forum,
I have been trying to solve a Power System Generation Expansion Optimization problem. Given that it is a loop that solves the mathematical problem, and then in a second phase executes some tests as a verifier, I built an optimizer object to repeatedly solve the similar problem. I want to constrain the addition of both solar and wind as to make additions in a not so disproportionate way, thus I want to limit that for any given year in the future, the total energy of any lump technology cannot be more than 3 times the amount of the other. I tried to resolve this by adding an implies with the following code:

d = binvar(2, study_horizon,'full'); % The study horizon is 30 years, so this variable is a matrix 2x30.
...
for t = 1:study_horizon      % the array VarConstraints is built only for the implies. When calling the optimizer command I call it like: optimizer(Constraints + VarConstraints,Cost,options,{...},{...})
   VarConstraints = VarConstraints + [sum(d(:,t)) == 1,    % This line says that the solver must choose one of the following cases
                               implies(d(1,t),[0 <= sum(Exist_git(solar,t)) + sum(Cand_git(solar,t)) <= 3*(sum(Exist_git(wind,t)) + sum(Cand_git(wind,t)))]);       % This lines can be read as: 0 <= sum(Energy from solar at year t) <= 3*sum(Energy from wind at year t)
                               implies(d(1,t),[0 <= sum(Exist_git(wind,t)) + sum(Cand_git(wind,t)) <= 3*(sum(Exist_git(solar,t)) + sum(Cand_git(solar,t)))])];      % This lines can be read as: 0 <= sum(Energy from wind at year t) <= 3*sum(Energy from solar at year t)
end

I am trying to add some bounds because I am seeing a Warning that reads: "Warning: You have unbounded variables in an implication leading to a lousy big-M relaxation. "

I am seeing that the decision variables involved are Exist_git and Cand_git, but I have a previous constraint over Exist_git and Cand_git that read:

for t = 1:study_horizon
    Constraints = Constraints + [(0 <= Exist_git(:,t) <= Exist_Upper_Limit): 'Existing Generator Bounds'];  % There is no decision regarding existing infrastructure
    Constraints = Constraints + [(0 <=  Cand_git(:,t) <=xit(:,t).*Candidate_Upper_Limit): 'Candidate Generator Bounds']; 
end

Xit is the integer decision variable about constructing or not technology i at year t. In case it decides to invest in that technology, it has the option to make integer amount of investments across time.
Now, we observe that one is explicitly bounded, while the other is only implicitly bounded. Nonetheless, even when we see the warning we are observing a good or plausible answer as reads:

As we can see in the image the quotient of Solar-Cand and Wind-Cand is about 3:1. So the solver decides to invests in 3 times many solar generators as wind.

However, to eliminate the Warning I add an explicit (trivial) bound over the decision variable Cand_git that reads:

for t = 1:study_horizon
     VarConstraints = VarConstraints + [0 <= Cand_git(:,t) <= Max_Install]; % 
end 

MaxInstall is an array that dictates the maximum amount we can install for a given technology. For instance, solar and wind have a top maximum total amount of 40,000 MW.
Moreover, I have some other constraints relating to this topic that read:

for t = 1:study_horizon
    Constraints = Constraints + [xit(:,t).*Cand_Upper <= Max_Install] % xit is the integer decision variable decision v
end

To be clear, again, xit is the integer decision variable to invest or not in any given technology. So, Cand_Upper is the Installed Capacity of a single unit of technology i. For example, I am modelling solar with an individual upper limit of 2,250 MW and with a Max Install total of 40,000MW. Other technologies have different parameters.

Coming back to the original point, when I add the trivial constraint over Cand_git, the implies warning disappears, but the solution is unexpected as shown in the following figure:

As you can see in year 2029, as well, in previous years, the quotient of Solar-Cand and Wind-Cand is about 7:1. Moreover, this is a very unexepected result as the added constraint was nothing but a triviality implied within other constraints. But it appears to not be respecting the implies constraints (VarConstraints)

Is there something I am misinterpreting? Is there something flawed within my code? I can share the code, but it is astoundingly long, although easy to follow as it is heavily commented.

Finally, I would like to ask what to do to be able to participate in the Google Forums group, as I initially intended to make this query there, but I learnt that I cannot open a new discussion there.

[johanlofberg]

Don't remove the explicit bound on Cand_git, it just leads to numerically bad models

You never use d(2,t) which is weird

[johanlofberg]

Also, when supplying examples to debug, you must remove everything not necessary for the bug to occur. As you only talk about wind and solar, remove everything else (while still seeing the thing you don't expect). Also, if it happens on one year, why include other years. For an outsider, this code is simply to complex to understand, you should be able to distill it significantly down to a minimal reproducible example. Remove simplify simplify

[johanlofberg]

If your issues is that the "Solar-cand" row should not be larger than 3x"wind-cand", it is not due to any bug in implies, as manual modelling yields the same

        aa = sum(Exist_git(Flags.flag_EcoThermal_exist+1:Flags.flag_EcoThermal_exist+first_wind_exist-1,:,t)*(8760*GEP_Load.LD)) + sum(Cand_git(Flags.flag_EcoThermal_cand+1:Flags.flag_EcoThermal_cand+first_wind_cand-1,:,t)*(8760*GEP_Load.LD)) ;
        bb = 3*(sum(Exist_git(Flags.flag_EcoThermal_exist+first_wind_exist:Flags.flag_renewables_exist,:,t)*(8760*GEP_Load.LD)) + sum(Cand_git(Flags.flag_EcoThermal_cand+first_wind_cand:Flags.flag_renewables_cand,:,t)*(8760*GEP_Load.LD)));
        cc = sum(Exist_git(Flags.flag_EcoThermal_exist+first_wind_exist:Flags.flag_renewables_exist,:,t)*(8760*GEP_Load.LD)) + sum(Cand_git(Flags.flag_EcoThermal_cand+first_wind_cand:Flags.flag_renewables_cand,:,t)*(8760*GEP_Load.LD));
        dd = 3*(sum(Exist_git(Flags.flag_EcoThermal_exist+1:Flags.flag_EcoThermal_exist+first_wind_exist-1,:,t)*(8760*GEP_Load.LD)) + sum(Cand_git(Flags.flag_EcoThermal_cand+1:Flags.flag_EcoThermal_cand+first_wind_cand-1,:,t)*(8760*GEP_Load.LD)));
        VarConstraints = VarConstraints +  [sum(d(:,t)) == 1;       
                    %    implies(d(1,t),[0 <= aa]);
                     %   implies(d(1,t),[aa <= bb]);
                      %%  implies(d(2,t),[0 <= cc])
                        implies(d(2,t),[cc <= dd])
                     aa >= -1e12*(1-d(1,t));
                     bb - aa >= -1e12*(1-d(1,t));
                     cc >= -1e12*(1-d(2,t));
                     dd - cc >= -1e12*(1-d(2,t));
                        ];

[mrkento11]

Thank you for your continuous help Professor.

I will go back to checking the mathematical formulation to check my mistake, and if I indeed translated correctly to Yalmip.

Moreover, I apologize for not simplifying enough. Will try to do a better job at simplification in case of future questions.

Thanks.

[johanlofberg]

Easiest is to skip optimizer construction and make call directly using optimize with added constraints on parameters, as you then directly can see if any constraint is violated (i.e. is there a bug in solver), and check the values of variables and constraint that you think are wrong. Using optimizer is hard to debug as it doesn't assign any values so you cannot check things easily

[mrkento11]

Just to let you know Professor, I found out the mistake and it was a reasoning mistake rather than a coding or software mistake (Yalmip).
I was experiencing unexpected results because as you point out I was seeing that "Solar-cand" should not be larger than 3*"wind-cand". I wanted to solar not grow disproportionately bigger than wind, so I tried to add a constraint for that. In the small test case I am using solar is much more attractive than wind, but I wanted to also consider if wind might be more attractive than solar. I thought that a good approach was to add a constraint for solar <= 3xWind, but also add in case wind is more attractive to be wind <= 3xSolar. {OR(solar <= 3xWind, wind <= 3xSolar)}.

I was using implies somewhere else within the whole code and thought that I should use it there as well. But, as you pointed out the manual modeling gave out the same answer. I observed that implies (obviously) decides to abide by the more lenient constraint among the two, and the constraint {wind <= 3*solar} was always being abided. Therefore, I realized that it was a reasoning mistake on my behalf to think that I need an implies constraint there (no implies needed), when I should add both constraints, and let the solver find out which of them both is the restricting constraint {AND(solar <= 3xWind, wind <= 3xSolar)}.

Thank you for your help on this topic, Professor.