I have a problem regarding modelling an objective function which should only count if it is positive. I am using a MIP model with binary variables. So I have the following equation (not in GAMS yet):
surplusPower(t) =PVgeneration(t) - electricalPowerTotal(t);
if (surplusPower(t) < 0) {
surplusPower(t) = 0;
}
The variable electricalPowerTotal(t) also depends on binary decision variables to schedule electrical loads, whereas PVgeneratrion (t) is predefined and fixed.
The goal is to minimze the sum (for all t) of surplusPower(t).
Does anyone have an idea how I could do that? I am struggeling with the if-condition. I’d appreciate your help.
This is a branch of optimization called “disjunctive programming”. Much can be said about it, but one of the easiest way to implement this for your case is to use an auxiliary binary variable in what is called a “big-M” approach. I’ve put together a small example for you, in this case a general approach that can be improved depending on the objective function. But, perhaps it will be enough for you to get the idea.
Thanks cladelpino for suggesting this interesting approach,
I have a question regarding it:
I do not really understand how the values for the variable isPositive(t) are determined? I understand the purpose of this binary variable but do I have to specify the values before the optimization (this would not be what I want) or is it a decision variable (meaning that the solver decides about its value)?
They are, as you call them, “decision variables”. Check out my ready to run example in the prior post, I’ve never set their values.
In that example, by using constraints bm1a,bm1b,bm2a,bm2b you can see that the feasible region of the problem is limited to having isPositive = 1 only when surplusPower > 0. The usual approach to study these kind of formulations is to use a “truth table”, where you write each constraint for values of isPositive = 0 and isPositive = 1 and look at them. You will instantly notice how they work.
Thanks a lot cladelpino for your answer and your help,
unfortunately I do not understand everything:
In that example, by using constraints bm1a,bm1b,bm2a,bm2b you can see that the feasible region of the problem is limited to having isPositive = 1 only when surplusPower > 0.
I see that if isPositive=1 → positiveSurplusPower(t) = surplusPower(t)
and if isPositive=0 → positiveSurplusPower(t) = 0
But I do not understand why isPositive = 1 when surplusPower > 0. Let’s say surplusPower(t) = -10 and isPositive(t) = 1:
bm1a(t): positiveSurplusPower(t) =l= -10
bm1b(t): positiveSurplusPower(t) =g= -10
→ positiveSurplusPower(t) = -10
bm2a(t): positiveSurplusPower(t) =l= bigM
bm2b(t): positiveSurplusPower(t) =g= -bigM
So these equations are satifsfied if we for instance set bigM=1000, because -10 =l=1000 and -10 =g= -1000.
Furthermore I did not really get your idea of the truth table (maybe this is related to my question above):
The usual approach to study these kind of formulations is to use a “truth table”, where you write each constraint for values of isPositive = 0 and isPositive = 1 and look at them. You will instantly notice how they work.
Would you mind saying some further words about that approach (I know what a truth table is, however I do not understand how to derive the constraints out of it).
But I do not understand why isPositive = 1 when surplusPower > 0. Let’s say surplusPower(t) = -10 and isPositive(t) = 1:
bm1a(t): positiveSurplusPower(t) =l= -10
bm1b(t): positiveSurplusPower(t) =g= -10
→ positiveSurplusPower(t) = -10
bm2a(t): positiveSurplusPower(t) =l= bigM
bm2b(t): positiveSurplusPower(t) =g= -bigM
This combination is infeasible because I declared positiveSurplusPower as a “positive variable” (line 4), which implicitly sets a lower bound of 0 for its value.
Here is a truth table for you (but I think that you’ve got it right now that we’ve clarified the implicit bound)
The following fact is important because it makes modeling easier:
In your GAMS model you could write:
positive variable surplusPower(t);
variable z objective variable;
[...]
equation eq(t)
obj objective function;
eq(t).. surplusPower(t) =g= PVgeneration(t) - electricalPowerTotal(t);
obj.. z =e= sum(t, surplusPower(t)) + [...];
[...]
model m /all/;
solve m use mip min z;
The minimization pushes surplusPower(t) down. Since it is defined as positive variable, it has a lower bound of 0.
Apparently, this only works if you are minimizing sum(t, surplusPower(t)).
first I want to try it with the bigM approach suggested by cladelpino (afterwards I will also try the other approach). So I inserted the following code:
The model is immediately solved and the objective value is 0. When you have a closer look at the results, you’ll see that the variable surplusPower(t) has the expected values (positive and negative values, depending on pvGenerationTotal(t) - electricalPowerTotal(t)) while the variable positiveSurplusPower(t) is always 0. So obviously there is a semantic mistake in the code, that I do not see. Does anyone have an idea?
( Disclaimer : Fred suggestion comes into the domain of what I meant with
(BigM is) a general approach that can be improved depending on the objective function
. I’m so used to not being able to make assumptions about the objective function that I remember bigMs better than I can remember the simple “linearizations”. If you don’t see an objective function change in the short term, using the “linearization” suggested by Fred is the way to go. Anyway, disjunctive models are a lot of fun, and bigM is a nice tool to learn and have around. )
Regarding your attempt at implementation: I don’t see any obvious error in your equations.
A common matter with these kind of formulations: What are the values of isPositive ? If you see small non zero values (on the order of 1/bigM), you may need to tighten the integer tolerance of your solver (option epint in cplex, I normally set it at “0” with no issues I could attribute to that so far)
the values for isPositive are 0 for every timeslot in the results.
There is one further aspect that I generally don’t understand: What is the difference (in the result file) between the values of a variable and a corresponding equation that defines this variable? Let’s have a look at the following variables that are defined via equations:
Following points are incomprehensible for me in the results:
the values for eq_pvGenerationTotal and pvGenerationTotal(t) are equal to the ones of pvGenerationTotal(t)
the variable surplusPower(t) has different values (positive and negative) depending on pvGenerationTotal(t) - electricalPowerTotal(t), but the equation eq_surplusPower(t) is always 0
the variable electricalPowerTotal(t) and the equation eq_electricalPowerTotal(t) are almost identical with the difference that in the variable some timeslots have an increased power values of 3000. Whenever a decision variable x(t, household) or y(t,household) has the value 1 (binary variables), 3000 is added which is the constant value for electricalPower(household) .
I don’t understand why there are differences among those different variables and equations. Maybe this has something to do with the problem that my model is not solves in a reasonible way (but this does not necessarily have to be the case).
It seems that the example from post #2 misses to properly link variables SurplusPower and isPositive. The example doesn’t work if the objective is minimized instead of maximized.
I inserted the two equations. Now a reasonable solution is found But only for one household. If I increase the number of households I always get the error message:
“Row ‘eq_electricalPowerTotal(144)’ infeasible, all entries at implied bounds”
The problematic equation looks like this (as posted before):
I do not understand this as I did not set any bound for the eq_electricalPowerTotal. Besides I am questioning, why this problem does not occur when only using 1 household?
By the way: I read the section on Solution Listing in the documentation but still I don’t understand how the values of equations in the solution are calculated. There it is only wirtten:
“Equation level .l Level of the equation in the current solution, equal to the level of all terms involving variables.”
It’s probably easier to help if you could post your entire code such that it can be executed.
As mentioned before, there are various formulations to retrieve only the positive part of a variable. One that appears a bit mor intuitive to me is
It seems that the example from post #2 misses to properly link variables SurplusPower and isPositive. The example doesn’t work if the objective is minimized instead of maximized.
Yes, my bad, sorry . This last formulation that Fred presents is much better.
Peter:
I do not understand this as I did not set any bound for the eq_electricalPowerTotal.
The message from the solver is telling you that the equation is infeasible and all of the values of variables appearing in it cannot be changed to make it feasible (because they have reached either a explicit or implicit bound), it is not making reference only to the equation “bound”.
On equation levels: Picture your equation in a way where every variable is on the left hand side of the operator (=g=,=l=, or =e=) , and on the right hand side you have (if there is any) the “constant” (determined only by parameters and constants) term. The equation level is the current value of the left hand side, considering the current values of the variables.
By the message that you posted I think that you are using CPLEX. The drawback is that you won’t get back the infeasible solution to ease your debug, but CPLEX has some conflict resolution capabilities. https://www.gams.com/latest/docs/S_CPLEX.html#CPLEXiis
@cladelpino
Thanks cladelpino for your answer and help,
On equation levels: Picture your equation in a way where every variable is on the left hand side of the operator (=g=,=l=, or =e=) , and on the right hand side you have (if there is any) the “constant” (determined only by parameters and constants) term. The equation level is the current value of the left hand side, considering the current values of the variables.
Unfortunately this is not possible as I am defining variables in the equations. So my equations are not only constraints but sometimes only define formulars for variables. For example:
By the message that you posted I think that you are using CPLEX. The drawback is that you won’t get back the infeasible solution to ease your debug, but CPLEX has some conflict resolution capabilities. > https://www.gams.com/latest/docs/S_CPLEX.html#CPLEXiis
I included the iis in the option file but I cannot understand anything from the additional output which looks like this
I used this approach and for one household everything is okay. But the same problem which occurs when using the other approaches with more than one houshold also occurs with this new appraoch. I get the message that the problem is integer infeasible.
We have discussed several topics in this thread and should try to keep things separate.
There was the original question how to consider only the positive part of a variable. Some bigM formulations have been discussed. If you understand the reformulation(s), I think this question has been answered.
There was the question on how to interpret equation levels in the solution listing.
Claudio’s explanation is perfectly fine:
I do not understand your objection. If an equation defines a variable (e.g. as sum of other variables) and there are no other constants involved, all variables are moved to the left hand side and the level of the equation is zero (in a feasible solution).
Your model is integer infeasible. It seems that you suspect the cause of infeasibility in an erroneous bigM reformulation. But this is not necessarily the cause. There can be many reasons why your model is infeasible. Unfortunately, you still have not shared your code to allow other users to reproduce this. That makes it difficult to provide target-oriented help for your particular model.
Your model is (integer) infeasible, so there is no feasible point to all the constraints and integrality conditions you have specified. That means, there is either a mistake implementing some constraint, or the combination of model and data really allows for no (integer) feasible solution. This state of the optimization model is in general not a bug or problem. It is one of the possible theoretic outcomes of an optimization run.
The real work it to figure out where the infeasibility comes from. There is no single way of doing this. It very much depends on your knowledge and understanding of the original problem and it’s implementation as a model. If you are convinced, that there must be a feasible solution to your problem, maybe you can provide such a solution by hand. If that is possible, a good way to go about analyzing infeasibilities is to provide such a “feasible” solution to the problem by manually setting the variable level values (x.l(i,j) = …) and then generating the model with a full equation listing (option limrow=1e9;) This will flag the constraints that are infeasible with your “feasible” solution.
If you don’t want to work with your own feasible solution, you can also ask Cplex you give you the smallest relaxation of your model to make it feasible (look for option FeasOpt: https://www.gams.com/latest/docs/S_CPLEX.html#CPLEXfeasopt ). This is just an alternative to the iis option mentioned by Claudio. What you posted is the log of the conflict refiner. It shows that Cplex found an irreducible infeasible set of 5 constraints/variable bounds. In the lst file these equations/variables should be reported between solve summary and solution listing.
I run into this thread and i found it useful but it unfortunately it doesn’t completely solve my problem.
I’m running an optimization problem where I have to minimize my objective function.
my objective function is, in a simplified version, the following one:
z =e = sum(t, sum(i,production_costs(i)*power_production(i,t)))
Now, I have to add an element to this that should look like this:
z =e =[...] - sum(j, TWC(j)*HS(j))
My problem is the following. TWC(j) (where j is a set of 5 different units) should be defined as follow (not in GAMS language):
What Fred told you is very correct. Just to sum up, would you please follow his recommendation and:
Provide a feasible solution for the problematic situation (more than one household) and use the equation listing to tell which constraints are infeasible.
At least I would expect you to do 1 and 2 on your own before I run your actual model.
Side note:
Unfortunately this is not possible as I am defining variables in the equations. So my equations are not only constraints but sometimes only define formulars for variables. For example:
This is a very common misconception. The sooner you forget about it, the better your path in optimization will be. There is no difference between one type of equality constraint or the other, neither is the problem solved sequentially.
The best model to reason about optimization problems is: Constraints are scalar functions involving variables. Their value is then restricted by an operator and a constant term (not involving ANY variable).
In the example that you posted, you are expressing that, in feasible solutions, only combinations of values that make the difference between
For this equation, the equation level is exactly the value of this difference, so if you see an constraint value of 1, this means that the constraint is currently infeasible, and in the current solution, electricPowerTotal(t) is greater by one unit than the summation.
You can use pretty much the same reformulation as Peter. The only difference is, that you are interested whether a variable is <= or >= 1000 instead of zero. Please consider the following modified example.
Your approach with HS.l in the dollar condition of an equation does not work for several reasons.
If you use the .l attribute of a variable in a model the refers to the level of that variable before the model is solved. Just dropping the .l is also no solution. Tying the generation of a constraint to a variable whose level is determined while the corresponding model is solved, does not fit the general concept of a Mathematical Program. Constraints have to be known before the model is solved
This is a very common misconception. The sooner you forget about it, the better your path in optimization will be. There is no difference between one type of equality constraint or the other, neither is the problem solved sequentially.
The best model to reason about optimization problems is: Constraints are scalar functions involving variables. Their value is then restricted by an operator and a constant term (not involving ANY variable).
So you suggest that I should not define variables in equations but rather define them immediately? So instead of having for example
For this equation, the equation level is exactly the value of this difference, so if you see an constraint value of 1, this means that the constraint is currently infeasible, and in the current solution, electricPowerTotal(t) is greater by one unit than the summation.
I am sorry but I have to admit that I don’t understand what you are trying to tell me. What do you mean with seeing an constraint value of 1? If the problem is infeasible (according to the solver) I do not see any constraint value because the results are not created. When there is a solution the result file is created by GAMS and there you most often find values greater than one for the =e= equations.
. This last formulation that Fred presents is much better.