Optimize reservoir hedging rules

GAMS
1

Hey guys，

We try to optimize the reservoir water supply hedging rule using a long series of runoff data. The rules take the form of reservoir storage corresponding to 24 periods of the year.When the initial reservoir storage is smaller than the corresponding time rule volume, the water supply is the coefficient multiplied by the water demand. When the initial reservoir storage is larger than the reservoir storage of rule corresponding time period, the water supply is the water demand.

I can only write the judgment between the first year reservoir storage and the rule, and when the input data becomes a long series of data for many years, it will not be processed. For example, the reservoir storage storage of the 25th time period in the long series should be compared with the first time period in the rule. The code I tried is as follows, but it reports an error.

Set n nodes
/Res_SH,Res_LW,Res_YNH,
HHK,SLZ,BDHK,CHK,WZD,HYT,
SLZ_pad,BDHK_urruin, BDHK_pad,BDHK_oth,YNH_pad,YNH_oth,YNH_agr,DL_ur/
;

Set Y years
/1956*2017/;
Alias(Y,Y1);

Set M months
*/1*24/
/Jan, Feb, Mar, Apr_first, Apr_mid, Apr_last, May_first, May_mid, May_last,
Jun_first, Jun_mid, Jun_last, Jul_first, Jul_mid, Jul_last, Aug_first, Aug_mid, Aug_last,
Sep_first, Sep_mid, Sep_last, Oct, Nov, Dec/
;
Set t  calculation step
/1*1488/
;

Positive variables
D(t,n)         Diversions water
R_pad(t,n)     Release for pad user
R_ur(t,n)      Release for urban user
R_oth(t,n)     Release for other user
R_agr(t,n)     Release for agricaltire user
R_eco(t,n)     Release for ecological quantity
A(t,n)         Abandon water quantity
S(t,n)
Rule(M,n)
;

Binary variables
Yp1(t)
Yp2(t)
;

Equations
* Constraints--Water balance for each node
B_ReYNH
B_YNHpad1
B_YNHpad2
B_YNHpad3
B_YNHpad4
;

B_ReYNH(t)..
S(t,'Res_YNH') =e= beg_S('Res_YNH')$(ord(t) EQ 1)
        +  S(t-1,'Res_YNH')$(ord(t) GT 1)
        +  Inflow(t,'Res_YNH')
        -  Evaporation(t,'Res_YNH')
        -  Leaking(t,'Res_YNH')
        -  R_pad(t,'YNH_pad')
        -  R_oth(t,'YNH_oth')
        -  R_agr(t,'YNH_agr')
        -  R_eco(t,'Res_YNH')
        -  A(t,'Res_YNH');

B_YNHpad1(t)..
Yp1(t)+Yp2(t)=1;
B_YNHpad2(t)..
R_pad(t,'YNH_pad')=e= Yp1(t)*Demand(t,'YNH_pad')*(1-Dam_dep('YNH_pad'))+Yp2(t)*Demand(t,'YNH_pad');
B_YNHpad3(t)..
S(t,'Res_YNH')=l= Rule(M,'YNH_pad')$(ord(t)=ord(M)+(ord(Y)-1)*24)*Yp1(t)+Vupper(t,'YNH_pad')*Yp2(t);
B_YNHpad4(t)..
S(t,'Res_YNH')=g= Vlower(t,'YNH_pad') * Yp1(t)+(Rule(M,'YNH_pad')$(ord(t)=ord(M)+(ord(Y)-1)*24)+10e-6)*Yp2(t);

where beg_S(‘Res_YNH’), Inflow(t,‘Res_YNH’),Evaporation(t,‘Res_YNH’),Leaking(t,‘Res_YNH’),Demand(t,‘YNH_pad’),Dam_dep(‘YNH_pad’) is the parameter
and Rule(M,n) is the decision variables we want to optimize.

The error is : Set for ‘ord’ is not controlled.

Thanks

2

The code you provide does not compile (because of other errors than you describe). Please provide a complete example that leads to your described problem.

GAMS is not so intuitive when it comes to sorting since it maintains (due to efficiency) a single list of labels. If set membership goes against this natural “label order” the set is considered unordered and you can’t do lag and leads or “ord”. There are usually ways around this (if you really need the set ordered). But to give you better advice, one needs to reproduce your issue. The sets in the model snippet you provide are all ordered.

-Michael

3

:smiley:

4

The complete code was uploaded with data and model descriptions.

Thanks again for you kind help！

5

There are plenty of compilation errors in the code. I interpreted the error message you posted earlier incorrectly, I thought this was about order, but this was about an uncontrolled index used in ord. You seem to struggle with the concept of “parallel” statements in GAMS. You control indexes by using them on the left (a(i,j,k) = … controls/iterates over sets/indexes i, j, and k). On the right you can have more indexes but they also need to be controlled, e.g. by a sum operator (a(i,j,k) = sum(l, …)). In this example you can use controlled indexes i,j,k, and l inside the sum but only i,j,k outside the sum. So looking at your first compilation error:

B_YNHpad7(t)$(ord(t)>1)..
S(t,'Res_YNH')=l= Rule(M,'YNH_pad')$(ord(M)=mod(ord(t),24)+1)*Yp1(t)+Vupper(t,'Res_YNH')*Yp2(t);

it is clear that you have difficulties with this concept. In your code the index M is not controlled. The “left hand side”, the equation definition controls index t but you just use index M on the right. What M should GAMS take? The $() condition seems to indicate that you want a particular M depending on t. The way to handle this is controlling M and then imposing the condition. Here is the sum is not really a sum but selects probably a single (if at all) M:

B_YNHpad7(t)$(ord(t)>1)..
S(t,'Res_YNH')=l= sum(M$(ord(M)=mod(ord(t),24)+1), Rule(M,'YNH_pad')*Yp1(t))+Vupper(t,'Res_YNH')*Yp2(t);

one could also write it like this

B_YNHpad7(t,M)$(ord(t)>1 and ord(M)=mod(ord(t),24)+1))..
S(t,'Res_YNH')=l= Rule(M,'YNH_pad')*Yp1(t)+Vupper(t,'Res_YNH')*Yp2(t);

If there is always (for ord(t)>1) a single M that matched the condition ord(M)=mod(ord(t),24)+1) then these two formulations are equivalent. If there is a chance that such an M does not exist for a given t then the prior solution still has the constraint “S(t,‘Res_YNH’) =l= Vupper(t,‘Res_YNH’)*Yp2(t);” while in the latter solution there is no constraint for such t.

Hope this helps,
-Michael

6

Hi Micheal:
Thank you very much for your answer.
Now the model runs, but it is infeasible.

The first thing I would like to inquire is what is the way to convert the nonlinear equation to linear, such as equation B_YNHpad3, B_YNHpad4, B_YNHpad7, B_YNHpad8.
Another thing I would like to inquire is how to output the set of infeasible equations and what causes the model to be infeasible?
I tried to use Baron as solver, but there’s no output for an IIS on my model.

I have attached the GAMS file and related excel file that requires to run the code. I will appreciate any help.

7

Your MINLP is actually a MIQCP and hence you can use Gurobi. With Gurobi’s IIS you get:

               S O L V E      S U M M A R Y

     MODEL   YNHWZD              OBJECTIVE  obj
     TYPE    MIQCP               DIRECTION  MINIMIZE
     SOLVER  GUROBI              FROM LINE  592

**** SOLVER STATUS     1 Normal Completion
**** MODEL STATUS      19 Infeasible - No Solution
**** OBJECTIVE VALUE               NA

 RESOURCE USAGE, LIMIT         25.969 10000000000.000
 ITERATION COUNT, LIMIT         0    2147483647
 EVALUATION ERRORS             NA             0
Gurobi full license.
Gurobi library version 9.5.1
Reading parameter(s) from "C:\Users\mbuss\Downloads\gurobi.opt"
>>  iis 1
Finished reading from "C:\Users\mbuss\Downloads\gurobi.opt"
--- GMO setup time: 3.41
Space for names approximately 4.20 Mb
MIQCP status(3): Model was proven to be infeasible.
Computing IIS...
An IIS is a set equations and variables (ie a submodel) which is
infeasible but becomes feasible if any one equation or variable bound
is dropped.

A problem may contain several independent IISs but only one will be
found per run.
Number of equations in the IIS: 2
B_YNHpad1 = 1
B_YNHpad3 < -22579.3
Number of variables in the IIS: 1
Rule(1,YNH_pad) < 22579.3

-Michael

8

Can you explain why the model is miqcp? There is no quadratic term in the model.
Also, when solving miqcp with cplex solver, the error is reported as follows.

CPLEX Error 5002: ‘QCP_row_for_B_YNHpad3’ is not convex.

The latest code and input files are attached
Thank you very much for your help!

9

Silence,

When I run your model with GAMS 39.2, I get this message in the log:

CPLEX Error 5002: ‘B_YNHpad3’ is not convex.

Looking at your GAMS source, I see equation B_YNHpad3:
B_YNHpad3… beg_S(‘Res_YNH’)=l= Rule(‘1’,‘YNH_pad’)*Yp1(‘1’)+Vupper(‘1’,‘Res_YNH’)*Yp2(‘1’);

The multiplication of variables Rule and Yp1 make this constraint quadratic.

-Steve

10

Hi Steve，
The error you mentioned I have noticed. So, we would like to get help on how to fix this error.Thank you