Maximizing the difference in the order of binary variables

PeterBe · August 2, 2017, 9:22am

Hi guys,

the problem I am facing is going to be difficult to explain and I strongly question whether something like this is solvable with GAMS. But I’ll try:

Lets say we have two binary decision variables x(t) and y(t) and lets assume they have the following values:
t: 0 1 2 3 4
x: 0 0 1 0 1
y: 1 0 0 1 0

So now I want to calculate the difference in the order of the value 1 between x(t) and y(t), meaning: the first time y(t)=1 is for t=0 and the first time x(t)=1 is for t=2. So the difference in time is 2-0= 2. The second time y(t)=1 is for t=3 and the second time x(t) = 1 is for t=4. So the difference in time is 4-3=1. The sum of differences is 3.

Now my question is, if it is possible to maximize this difference in time by using GAMS with a solver. This problem also has some constraints that have to be taken into account.

Manassaldi · August 2, 2017, 3:43pm

Hi, I think you should use an auxiliary variable to define the first time of X and Y (same to the second time and third time and so on)

I think this can work :

eqfirstx(t).. 1 - x(t) + sum(tp$(ord(tp) lt ord(t)),x(tp)) + xfirst(t) =g= 1;
eqfirsty(t).. 1 - y(t) + sum(tp$(ord(tp) lt ord(t)),y(tp)) + yfirst(t) =g= 1;
eqdiffirst..                                   firstdif =e= sum(t,xfirst(t)*(ord(t)-1)) - sum(t,yfirst(t)*(ord(t)-1));
*Similarly for the second time of x:
eqsecondx(t,tp)$(ord(tp) lt ord(t)).. 1 - x(t) + 1 - x(tp) + sum(tpp$(ord(tpp) lt ord(t) and ord(tpp) ne ord(tp)),x(tpp)) + xsecond(t) =g= 1;
eqsecondy(t,tp)$(ord(tp) lt ord(t)).. 1 - y(t) + 1 - y(tp) + sum(tpp$(ord(tpp) lt ord(t) and ord(tpp) ne ord(tp)),y(tpp)) + ysecond(t) =g= 1;
eqdifsecond..                                   seconddif =e= sum(t,xsecond(t)*(ord(t)-1)) - sum(t,ysecond(t)*(ord(t)-1));
eqobj.. z =e=  firstdif + seconddif;

In “sum(t,xfirst(t)*(ord(t)-1)” I use “ord(t)-1” because “ord(t)” indicates the position of the set not its value

Manassaldi · August 2, 2017, 3:51pm

I forgot this equations

uniquefirstx..     sum(t,xfirst(t)) =e= 1;
uniquefirsty..     sum(t,yfirst(t)) =e= 1;
uniquesecondx.. sum(t,xsecond(t)) =e= 1;
uniquesecondy.. sum(t,ysecond(t)) =e= 1;

Bye!

PeterBe · August 3, 2017, 4:10pm

Hi Manassaldi,

thank you for your answer. Unfortunately I’m having a hard time understading what you posted. Even the first equation is not understandable for me (I have not tried to understand the other equations yet, because I am struggeling even with the first one):

eqfirstx(t).. 1 - x(t) + sum(tp$(ord(tp) lt ord(t)),x(tp)) + xfirst(t) =g= 1;

Let’s have a look a the exampe:
t: 0 1 2 3 4
x: 0 0 1 0 1
y: 1 0 0 1 0

For t=2
1- 1 + 0 + xfirst(2) =1 → xfirst(2) = 1

So xfirst(2) should be one. But this cannot be true, since xfirst has to have the value 2, according to my understandings.

Manassaldi · August 3, 2017, 7:24pm

Hi, “xfirst” identifies the position of the first time of x and the term “sum(t,xfirst(t)*(ord(t)-1))” corresponds to its order.

For your example:
t: 0 1 2 3 4
x: 0 0 1 0 1
y: 1 0 0 1 0

For t=2 (correspond to position 3 of the set)
eqfirstx(t=2) → 1- 1 + 0 + xfirst(2) >=1 → xfirst(2) = 1
For t=0 (correspond to position 1 of the set)
eqfirsty(t=0) → 1 - 1 + 0 + yfirst(0) >=1 → yfirst(0) = 1
If you expand "sum(t,xfirst(t)(ord(t)-1))" → xfirst(0)(1-1) + xfirst(1)(2-1) + xfirst(2)(3-1) + xfirst(3)(4-1) + xfirst(4)(5-1)
but only xfirst(2) take a value of 1, so: sum(t,xfirst(t)(ord(t)-1)) = xfirst(2)(3-1) = 1*(3-1) = 2
Similarly for y: sum(t,yfirst(t)(ord(t)-1)) = xfirst(0)(1-1) = = 1*(1-1) = 0
finally
eqdiffirst → firstdif = 2 - 0 = 2;

firstdif is the difference in the order of the value 1 between x(t) and y(t) for the first time

PeterBe · August 4, 2017, 12:34pm

Thanks once again Manassaldi for your good answer,

Now I understood how your concept works. You said that I’ll have to do this for every occurence of x(t)=1 and y(t)=1 I have to comments/questions regarding this:

The number of timeslots, inwhich x(t) = 1and y(t) =1 is not known in advance and has to be determined by the solver itself. As I said earlier there are some constraints that have to be taken into account.
Is there a way how to do this more “conveniently” considering that I’ll have to do this for 1440 timeslots. Or do I really have to use thousands of equations for that?

Manassaldi · August 4, 2017, 5:10pm

Hi, analyzing more in detail the problem maybe there is a better way to express it.
Perhaps there is a compact way of expressing the constraints.

I think you must to propose the maximum time in which x(t) = 1 and y(t) =1 before, but you have to modifies some restriction.
To consider the case where a “time slot” does not exist all these restrictions have to change to “=l=”
uniquefirstx… sum(t,xfirst(t)) =e= 1; → uniquefirstx… sum(t,xfirst(t)) =l= 1;
uniquefirsty… sum(t,yfirst(t)) =e= 1; → uniquefirsty… sum(t,yfirst(t)) =e= 1;
uniquesecondx… sum(t,xsecond(t)) =e= 1; → uniquesecondx… sum(t,xsecond(t)) =e= 1;
uniquesecondy… sum(t,ysecond(t)) =e= 1; → uniquesecondy… sum(t,ysecond(t)) =e= 1;

best regard

PeterBe · August 7, 2017, 8:08am

Thanks Manassaldi for your answer (and your general great effort ),

I have to admit that I do not really understand the things you said. What is meant by:

I think you must to propose the maximum time in which x(t) = 1 and y(t) =1 before, but you have to modifies some restriction.

As I said in the earlier post. The number of time in which x(t) = 1 and y(t) =1 is not known in advance. I could specify an upper bound for that, but the basic idea is that this number is not known. Are you assuming that there is no way how to handle that?

And you also said that I should change the constraints like this:

uniquefirstx… sum(t,xfirst(t)) =e= 1; → uniquefirstx… sum(t,xfirst(t)) =l= 1;
uniquefirsty… sum(t,yfirst(t)) =e= 1; → uniquefirsty… sum(t,yfirst(t)) =e= 1;

Why do I only have to change uniquefirstx and not uniquefirsty?

Manassaldi · August 7, 2017, 11:14am

Hi, sorry, it’s a mistake
uniquefirstx… sum(t,xfirst(t)) =e= 1; → uniquefirstx… sum(t,xfirst(t)) =l= 1;
uniquefirsty… sum(t,yfirst(t)) =e= 1; → uniquefirsty… sum(t,yfirst(t)) =l= 1;
uniquesecondx… sum(t,xsecond(t)) =e= 1; → uniquesecondx… sum(t,xsecond(t)) =l= 1;
uniquesecondy… sum(t,ysecond(t)) =e= 1; → uniquesecondy… sum(t,ysecond(t)) =l= 1;

Manassaldi · August 7, 2017, 1:50pm

Hi again, you could try with a BIGM reformulation for simplify the wole formulation of the problem
I think it’s better to start with 1 instead of 0
Try the next code:

set
t /1*1440/
;
alias(t,tp,tpp);
binary variable
x(t)
y(t)
;
variable
z
;
equation
eq1,eq2,eq3,eq4,eq5,eq6,eqobj
;

eq1(t,tp).. sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= ord(tp) + (1-xpos(tp,t))*1440;
eq2(t,tp).. sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= ord(tp) - (1-xpos(tp,t))*1440;
eq3(tp)..   sum(t,xpos(tp,t)) =l= 1;
*same for y
eq4(t,tp).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =l= ord(tp) + (1-ypos(tp,t))*1440;
eq5(t,tp).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =g= ord(tp) - (1-ypos(tp,t))*1440;
eq6(tp)..   sum(t,ypos(tp,t)) =l= 1;

eqobj.. z =e=  sum(tp, sum(t,xpos(tp,t)*ord(tp)) - sum(t,ypos(tp,t)*ord(tp)));

For example:
t: 1 2 3 4 5 6 … … 1440
x: 0 0 1 0 1 0 …(all 0) 0
y: 1 0 0 1 0 0 …(all 0) 0

now, only this variables take a value of 1:
xpos(‘1’,‘3’)=1 and xpos(‘2’,‘5’)=1
ypos(‘1’,‘1’)=1 and ypos(‘2’,‘4’)=1

you must check this but I think this can works

Best regard

PeterBe · August 9, 2017, 9:06am

Thanks Manassaldi for your answer,

but I do nor understand your code:

set
t /1*1440/
;
alias(t,tp,tpp);
binary variable
x(t)
y(t)
;
variable
z
;
equation
eq1,eq2,eq3,eq4,eq5,eq6,eqobj
;

eq1(t,tp).. sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= ord(tp) + (1-xpos(tp,t))*1440;
eq2(t,tp).. sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= ord(tp) - (1-xpos(tp,t))*1440;
eq3(tp)..   sum(t,xpos(tp,t)) =l= 1;
*same for y
eq4(t,tp).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =l= ord(tp) + (1-ypos(tp,t))*1440;
eq5(t,tp).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =g= ord(tp) - (1-ypos(tp,t))*1440;
eq6(tp)..   sum(t,ypos(tp,t)) =l= 1;

eqobj.. z =e=  sum(tp, sum(t,xpos(tp,t)*ord(tp)) - sum(t,ypos(tp,t)*ord(tp)));

I have several questions and comments:

xpos(tp,t) and ypos(tp,t) are not declared in your model description. They are only unsed in the equations. Is it not necessary to define them?
Basically I have to admit that I do not understand much of your code and I am confused (I spend many hours trying to understand). What are eq1, eq2 and eq3 ensuring? On the left side of eq1 and eq2 your are counting the number of times x(t)=1 until t. But on the right side you are using tp which I don’t understand. When xpos(tp,t) =0, eq1 and eq2 ensure that the sum equals ord(tp). Let’s say tp=1000 and t is any number (lets say t=1010) and xpos(tp,t)=1, then eq1 and eq2 - according to my understandig - force the sum of x(t)=1 for t in [0,…, 1005] to be 1000. That means that x(t) is almost always 1. But when I have constraints that forbid something like this there will not be a feasible solution.
I guess I have to insert equations for every timeslot and every binary variable, meaning that I have to insert - taking your suggestion - at least 4*1440 = 5760 equations manually. This is also not really doable for me. It would take a fairly long time and there will surely be mistakes in several of those equations which I’ll not be able to find.

I know that this question is a fairly difficult one. But I at least want to try to solve it. This is why I really appreciate your input (and the inputs of others)

Manassaldi · August 9, 2017, 2:09pm

Hi Peter,
I forgot to define the variables because I wrote them fast in the browser without checking them in GAMS.
xpos(tp,t) and ypos(tp,t) are binary variable

For example,
if xpos(‘1’,‘4’)=1 this mean thats the first time of x it happen in the fouth order (tp=1 and t=4)
if xpos(‘2’,‘10’)=1 this mean thats the second time of x it happen in the tenth order (tp=2 and t=10)
.
.
.
if xpos(‘50’,‘1000’)=1 this mean thats the fiftieth time of x it happen in the thousandth order (tp=50 and t=1000)

But, if for the rest of restriction x it happen only 50 time, the rest of the variables xpos take a values of 0
( xpos(‘51’,t)=xpos(‘52’,t)=xpos(‘53’,t)=…=xpos(‘1440’,t)=0 )

For the previusly example, xpos(‘1’,‘4’)=1 (tp=1 and t=4)
eq1(t,tp)… sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= ord(tp) + (1-xpos(tp,t))*1440;
eq2(t,tp)… sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= ord(tp) - (1-xpos(tp,t))*1440;
replacing and expanding
eq1(‘4’,‘1’)… x(‘1’) + x(‘2’) + x(‘3’) + x(‘4’) =l= 1 + (1-1)*1440;
eq2(‘4’,‘1’)… x(‘1’) + x(‘2’) + x(‘3’) + x(‘4’) =g= 1 - (1-1)*1440;
and
x(‘1’) + x(‘2’) + x(‘3’) + x(‘4’) = 1

if xpos(‘1’,‘4’)=1 the rest of xpos(‘1’,t) (t ne 4) take a value of 0, so
eq1(t,‘1’)… sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= 1 + (1-0)*1440;
eq2(t,‘1’)… sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= 1 - (1-0)*1440;
finally
eq1(t,‘1’)… sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= 1441;
eq2(t,‘1’)… sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= 1441;
eq1 and eq2 always be satisfied in this case

Analyzing your questions, I noted that you must to aggregate the eq7 and eq8
I think that this equations can represent all your problem. This focus is so-called “BigM relaxation”

Anyway, only you know all the problem and you must to vericated all equations.

I hope this can help you
bye!

set
t /1*1440/
;
alias(t,tp,tpp);
binary variable
x(t)
y(t)
xpos(tp,t)
ypos(tp,t)

;
variable
z
;
equation
eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eqobj
;

eq1(t,tp)..  sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= ord(tp) + (1-xpos(tp,t))*1440;
eq2(t,tp)..  sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= ord(tp) - (1-xpos(tp,t))*1440;
eq3(tp)..    sum(t,xpos(tp,t)) =l= 1;
*same for y
eq4(t,tp).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =l= ord(tp) + (1-ypos(tp,t))*1440;
eq5(t,tp).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =g= ord(tp) - (1-ypos(tp,t))*1440;
eq6(tp)..   sum(t,ypos(tp,t)) =l= 1;

eq7(tp,t)..  1-xpos(tp,t) + x(t) =g= 1;
eq8(tp,t)..  1-ypos(tp,t) + y(t) =g= 1;

eqobj.. z =e=  sum(tp, sum(t,xpos(tp,t)*ord(tp)) - sum(t,ypos(tp,t)*ord(tp)));

Manassaldi · August 10, 2017, 1:22pm

I think that the first formulation of the problem is better than the BigM.
In the first approximation you must to type more equation but the final size is lower.
Maybe, if you could found a compact way to express the first formulation the model will result small an easy to implement.

If an upper bound is proposed for the number of time in which x=1 (an upper bound not the exact number of time in which x=1) you can reduce the model size significantly.
For example: maximum number of time in which x=1 → 50

For the first formulation:

eqfirstx(t).. 1 - x(t) + sum(tp$(ord(tp) lt ord(t)),x(tp)) + xfirst(t) =g= 1;
eqsecondx(t,tp)$(ord(tp) lt ord(t)).. 1 - x(t) + 1 - x(tp) + sum(tpp$(ord(tpp) lt ord(t) and ord(tpp) ne ord(tp)),x(tpp)) + xsecond(t) =g= 1;
eqthridx(t,tp,tpp)$(ord(tp) lt ord(t) and ord(tpp) lt ord(t)).. 1 - x(t) + 1 - x(tp) + 1 - x(tpp) + sum(tppp$(ord(tppp) lt ord(t) and ord(tppp) ne ord(tp)  and ord(tppp) ne ord(tpp)),x(tppp)) + xthrid(t) =g= 1;
.
.
.
eqfiftiethx(t,tp)  --> Analogously to the second and thrid but taken from a 50

For the BigM formulation:

Parameter
Maxim /50/
;
eq1(t,tp)$(ord(tp) le Maxim)..  sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =l= ord(tp) + (1-xpos(tp,t))*1440;
eq2(t,tp)$(ord(tp) le Maxim)..  sum(tpp$(ord(tpp) le ord(t)),x(tpp)) =g= ord(tp) - (1-xpos(tp,t))*1440;
eq3(tp)$(ord(tp) le Maxim)..    sum(t,xpos(tp,t)) =l= 1;
*same for y
eq4(t,tp)$(ord(tp) le Maxim).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =l= ord(tp) + (1-ypos(tp,t))*1440;
eq5(t,tp)$(ord(tp) le Maxim).. sum(tpp$(ord(tpp) le ord(t)),y(tpp)) =g= ord(tp) - (1-ypos(tp,t))*1440;
eq6(tp)$(ord(tp) le Maxim)..   sum(t,ypos(tp,t)) =l= 1;

eq7(tp,t)$(ord(tp) le Maxim)..  1-xpos(tp,t) + x(t) =g= 1;
eq8(tp,t)$(ord(tp) le Maxim)..  1-ypos(tp,t) + y(t) =g= 1;

eqobj.. z =e=  sum(tp$(ord(tp) le Maxim), sum(t,xpos(tp,t)*ord(tp)) - sum(t,ypos(tp,t)*ord(tp)));

I hope this can help you
bye!

PeterBe · August 10, 2017, 3:07pm

Hi Manassaldi, Thank you for your answers

I still have questions:

In your first formulation I think you have to insert equations for all timeslots and the number of time in which x is has to be know in advance. Is that correct?

eqfirstx(t).. 1 - x(t) + sum(tp$(ord(tp) lt ord(t)),x(tp)) + xfirst(t) =g= 1;
eqfirsty(t).. 1 - y(t) + sum(tp$(ord(tp) lt ord(t)),y(tp)) + yfirst(t) =g= 1;
eqdiffirst..                                   firstdif =e= sum(t,xfirst(t)*(ord(t)-1)) - sum(t,yfirst(t)*(ord(t)-1));
*Similarly for the second time of x:
eqsecondx(t,tp)$(ord(tp) lt ord(t)).. 1 - x(t) + 1 - x(tp) + sum(tpp$(ord(tpp) lt ord(t) and ord(tpp) ne ord(tp)),x(tpp)) + xsecond(t) =g= 1;
eqsecondy(t,tp)$(ord(tp) lt ord(t)).. 1 - y(t) + 1 - y(tp) + sum(tpp$(ord(tpp) lt ord(t) and ord(tpp) ne ord(tp)),y(tpp)) + ysecond(t) =g= 1;
eqdifsecond..                                   seconddif =e= sum(t,xsecond(t)*(ord(t)-1)) - sum(t,ysecond(t)*(ord(t)-1));
eqobj.. z =e=  firstdif + seconddif;

Just a general question: x(t) and y(t) should still be binary decision variables. Ist it possible to addconstraints that only refer to x(t) and y(t) and not to xpos(tp, t) and ypos (tp, t)? For example the one I asked two weeks ago in this forum (https://newforum.gams.com/t/ensuring-realations-between-decision-variables/1905/1) and other constraints

Manassaldi · August 10, 2017, 5:14pm

Hi,

In your first formulation I think you have to insert equations for all timeslots and the number of time in which x is has to be know in
advance. Is that correct?

You have to propose the maximum number of time in which x=1, this is an upper bound. You don’t has to know the number of time in which x=1 in advance, just only an upper bound.
So, if this upper bound is equal to 50 you have to insert 50 equation for x and 50 for y.
Later, if x is equal 1 only forty time, xfortyone(t), xfortytwo(t),… and xfifty(t) will be equal zero.

2) Just a general question: x(t) and y(t) should still be binary decision variables. Ist it possible to addconstraints that only refer to x(t) and y(t) and not to xpos(tp, t) and ypos (tp, t)? For example the one I asked two weeks ago in this forum (viewtopic.php?f=9&t=10040) and other constraints

xpos(tp, t) and ypos(tp,t) are auxiliary variable only to identifies the order
Anyway, i suggest to you that use the first approximation

I hope this can help you
bye!

PeterBe · August 11, 2017, 8:05am

Hi Manassaldi,

So, if this upper bound is equal to 50 you have to insert 50 equation for x and 50 for y.

I don’t really have an upper bound. The only upper bound I have is the number of timeslots,which is 1440. So I have to insert 1440 *2= 2880 equations when using the first approximation? What about the BigM-relaxation? Do I also have to use thousands of equations there?

Manassaldi · August 11, 2017, 12:39pm

Hi Peter,
no, for bigM reformultation you don’t have to use thousands of equations. The model that I showed was complete, but when is expanded is very bigger.
As i say previously, if you could found a compact way to express the first formulation of the problem the model will result more small and easy to implement.

Bye!

PeterBe · August 14, 2017, 2:27pm

Hi Manassaldi,

i tried to execute the BigM-Relaxation approach of my model. However my computer was out of memory. The workspace was more than 27.000 Mb and then GAMS stopped with the message “Out of memory”.

I think by memory GAMS is referring to RAM and not to harddrive. Because on my harddrive there is more than 100 GB left but my RAM is only 8 GB.

Manassaldi · August 14, 2017, 3:25pm

Yes, this BigM formulation generates a very large model…

PeterBe · August 15, 2017, 7:41am

Yes, this BigM formulation generates a very large model…

Two question regarding this:

Does the first approximation generate a smaller model which does not need to much memory? Or do you know any other way how to reduce the required memory size?
Is there a way to roughly estimate how much memory this BigM formulation of my model needs? Because I could use other computers which have more than 8GB memory. But I should only use them if I know that it will be useful