XiongPengNUS / rsome

Robust Stochastic Optimization Made Easy
GNU General Public License v3.0
282 stars 54 forks source link

Numerical problems with equality constraints with rvar #46

Closed sometimesstudy closed 1 year ago

sometimesstudy commented 1 year ago

Hi,

When I actually applied RSOME to solve adaptive robust, I found that when there is an equality constraint on rvar, it is easy to be infeasible.

Take the adaptive robust power system scheduling as an example,

For shortly

  1. PV load is defined as PV = ldr(24)
  2. PV load fluctuation is uncertainty variable delta = rvar(24) PV.adapt(delta)
  3. uncertainty set is defined as zset = (delta>0 , delta<1, delta.sum()<=12)
  4. constraint PV = PV_pre +delta
    Pg + PV = PL

A deterministic model, or a adaptive robust model that does not contain equality constraints is feasible, but become infeasible when equality constraints exist.

Moreover, I have tried to rewrite the inequality constraint in the Adaptive Robust Lot-Sizing example into an equality constraint model.st(d <= y.sum(axis=0) - y.sum(axis=1) + x)->model.st(d == y.sum(axis=0) - y.sum(axis=1) + x)

It also becomes infeasible. So what can I do when the model contains equality constraints?

GalPerelman commented 1 year ago

@sometimesstudy I'm not sure if this is the cause of the infeasibility but it seems that the adaption is not as it should be As it is written now: PV.adapt(delta), all the decision variables are adapting all the random variables In time dynamic problems (such as economic dispatch) the LDR should be from a non-anticipative type This means that LDR for each decision variable contains only information from periods that are previous to the decision itself The Joint Production-Inventory example is very recommended for that

sometimesstudy commented 1 year ago

Hi, @GalPerelman

Thanks for your reply.

I modified PV.adapt(delta) to

for t in range(T_len):
    PV[t].adapt(delta[t])

But the model is still infeasible.

Indeed, PV is a LDR with length of 24, and delta is a RVAR with length of 24, denoting the fluctuation of Photovoltaic.

The relationship of PV and delta is

PV = PV0+delta

where PV0 is a given parameter.

chenzjames commented 1 year ago

In general, we shall avoid “equality” in a robust constraint, which could be problematic. Most cases, a weak inequality (either >= or <=) will do.

On 6 May 2023, at 15:32, Yuxuan @.***> wrote:

Hi, @GalPerelman https://github.com/GalPerelman Thanks for your reply.

I modified PV.adapt(delta) to

for t in range(T_len): PV[t].adapt(delta[t]) But the model is still infeasible.

Indeed, PV is a LDR with length of 24, and delta is a RVAR with length of 24, denoting the fluctuation of Photovoltaic.

The relationship of PV and delta is

PV = PV0+delta

where PV0 is a given parameter.

— Reply to this email directly, view it on GitHub https://github.com/XiongPengNUS/rsome/issues/46#issuecomment-1537078692, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWDYEQE2N5TKCKBW5MI3XEX5B3ANCNFSM6AAAAAAXX5A7JU. You are receiving this because you are subscribed to this thread.

GalPerelman commented 1 year ago

@sometimesstudy I think that it should be:

for t in range(1, T_len):
    PV[t].adapt(delta[:t])           # using the ":" before t means adaptation all until t

As I wrote above, I'm really not sure this is the cause of the infeasibility

sometimesstudy commented 1 year ago

@GalPerelman

Thanks for your help, but the model is still infeasible. And the following is my complete program.

Defined parameters


# import module
from gurobipy import *
import numpy as np
import pandas as pd
from rsome import ro                            # import the ro module
from rsome import grb_solver as grb 

# Defined parameters
dT = 1
T_set = 24
PG_max = 800  # DG出力上限
PG_min = 80  # DG出力下限
a = 0.67  # 成本系数
b = 0  # 成本系数
K_S = 0.38  # 单位充放电成本
PS_max = 500  # 充放电功率上限
ES_max = 1800  # 荷电状态上限
ES_min = 400  # 荷电状态下限
ES0 = 1000  # 初始荷电状态
eta = 0.95  # 储能充放电效率
K_DR = 0.32  # 单位调度成本
D_DR = 2940  # 总用电需求
D_DR_max = 200  # 最大用电需求
D_DR_min = 50  # 最小用电需求
P_DR0=[80,70,60,50,70,70,90,100,120,150,170,200,140,100,100,120,140,150,190,200,200,190,100,80]
PM_max = 1500  # 最大交互功率
lambda1 = [0.48, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48, 0.9, 1.35, 1.35, 1.35,
           0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 1.35, 1.35, 1.35, 1.35, 1.35, 0.48]
P_PV0=np.array([0,0,0,0,0,0.0465,0.1466,0.3135,0.4756,0.5213,0.6563,1,0.7422,0.6817,0.4972,0.4629,
       0.2808,0.0948,0.0109,0,0,0,0,0])*1500
P_L0=np.array([0.4658,0.4601,0.5574,0.5325,0.5744,0.6061,0.6106,0.6636,0.741,0.708,0.7598,0.8766,0.7646,
    0.7511,0.6721,0.5869,0.6159,0.6378,0.6142,0.6752,0.6397,0.5974,0.5432,0.4803])*1200

# Uncertainty
delta_u_PV=0.15
delta_u_L=0.1
tau_L = 12
tau_PV = 6
model = ro.Model('MGs')
Ppv_t = model.ldr(T_set)
B_pv_t = model.rvar(T_set)
for t in range(1,24):
    Ppv_t[t].adapt(B_pv_t[:t])
zset = (B_pv_t>=0, B_pv_t<=1, B_pv_t.sum()<=6)

# First stage variables
Us_t = model.dvar(T_set,vtype='B')
Um_t = model.dvar(T_set,vtype='B')

# Second stage variables
Pg_t = model.dvar(T_set)
Ps_ch_t = model.dvar(T_set)
Ps_dis_t = model.dvar(T_set)
Pdr_t = model.dvar(T_set)
Pdr1_t = model.dvar(T_set)
Pdr2_t = model.dvar(T_set)
Pbuy_t = model.dvar(T_set)
Psell_t = model.dvar(T_set)
Pl_t = model.dvar(T_set)
B_l_t = model.dvar(T_set)

# Constraints
obj1 = a*Pg_t+b
obj2 = K_S*(Ps_dis_t*(1/eta)+Ps_ch_t*eta)
model.st(Ps_dis_t <= Us_t*PS_max)
model.st(Ps_ch_t <= (1-Us_t) *PS_max)
model.st(eta*Ps_ch_t.sum()-1/eta*Ps_dis_t.sum() == 0)
model.st((ES0+eta*Ps_ch_t*dT-1/eta*Ps_dis_t*dT)>= ES_min)
model.st((ES0+eta*Ps_ch_t*dT-1/eta*Ps_dis_t*dT)<= ES_max)
obj3 = K_DR*(Pdr1_t+Pdr2_t)
model.st(Pdr_t.sum()==D_DR)
model.st(Pdr_t-np.array(P_DR0)+Pdr1_t-Pdr2_t==0)
obj4 = np.array(lambda1)*(Pbuy_t-Psell_t)
model.st(Pbuy_t-Psell_t==Ps_ch_t+Pdr_t+Pl_t-Pg_t-Ps_dis_t-Ppv_t)
model.st(Pbuy_t<=Um_t*PM_max)
model.st(Psell_t<=(1-Um_t)*PM_max)

# constraints of LDR and RVAR
model.st(Ppv_t==P_PV0+B_pv_t)
model.st(Pl_t==P_L0+B_l_t*delta_u_L)

# Objctive
model.minmax(obj1.sum()+obj2.sum()+obj3.sum()+obj4.sum(),zset)

# Solve
model.solve(grb)
print(model.get())
sometimesstudy commented 1 year ago

Hi, Dr. @chenzjames

I tried to convert the equality constraints into inequality constraints before, the model is feasible, but the optimization results are far from the results I calculated with the CCG algorithm. How can I avoid the gap?

In general, we shall avoid “equality” in a robust constraint, which could be problematic. Most cases, a weak inequality (either >= or <=) will do. On 6 May 2023, at 15:32, Yuxuan @.***> wrote: Hi, @GalPerelman https://github.com/GalPerelman Thanks for your reply. I modified PV.adapt(delta) to for t in range(T_len): PV[t].adapt(delta[t]) But the model is still infeasible. Indeed, PV is a LDR with length of 24, and delta is a RVAR with length of 24, denoting the fluctuation of Photovoltaic. The relationship of PV and delta is PV = PV0+delta where PV0 is a given parameter. — Reply to this email directly, view it on GitHub <#46 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWDYEQE2N5TKCKBW5MI3XEX5B3ANCNFSM6AAAAAAXX5A7JU. You are receiving this because you are subscribed to this thread.

chenzjames commented 1 year ago

Is it because that RSOME use LDR “approximation”, while CCG yields the exact solution?

On 6 May 2023, at 16:03, Yuxuan @.***> wrote:

Hi, Dr. @chenzjames https://github.com/chenzjames I tried to convert the equality constraints into inequality constraints before, the model is feasible, but the optimization results are far from the results I calculated with the CCG algorithm. How can I avoid the gap?

In general, we shall avoid “equality” in a robust constraint, which could be problematic. Most cases, a weak inequality (either >= or <=) will do. … <x-msg://6/#> On 6 May 2023, at 15:32, Yuxuan @.***> wrote: Hi, @GalPerelman https://github.com/GalPerelman https://github.com/GalPerelman Thanks for your reply. I modified PV.adapt(delta) to for t in range(T_len): PV[t].adapt(delta[t]) But the model is still infeasible. Indeed, PV is a LDR with length of 24, and delta is a RVAR with length of 24, denoting the fluctuation of Photovoltaic. The relationship of PV and delta is PV = PV0+delta where PV0 is a given parameter. — Reply to this email directly, view it on GitHub <#46 (comment) https://github.com/XiongPengNUS/rsome/issues/46#issuecomment-1537078692>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWDYEQE2N5TKCKBW5MI3XEX5B3ANCNFSM6AAAAAAXX5A7JU. You are receiving this because you are subscribed to this thread.

— Reply to this email directly, view it on GitHub https://github.com/XiongPengNUS/rsome/issues/46#issuecomment-1537084378, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWD7LQOFJQT4RZJ3YEILXEYAUZANCNFSM6AAAAAAXX5A7JU. You are receiving this because you were mentioned.

sometimesstudy commented 1 year ago

Dr. @chenzjames

Thanks a lot, and I wonder how to assign vtype='B' to rvar.

Is it because that RSOME use LDR “approximation”, while CCG yields the exact solution? On 6 May 2023, at 16:03, Yuxuan @.> wrote: Hi, Dr. @chenzjames https://github.com/chenzjames I tried to convert the equality constraints into inequality constraints before, the model is feasible, but the optimization results are far from the results I calculated with the CCG algorithm. How can I avoid the gap? In general, we shall avoid “equality” in a robust constraint, which could be problematic. Most cases, a weak inequality (either >= or <=) will do. … <x-msg://6/#> On 6 May 2023, at 15:32, Yuxuan @.> wrote: Hi, @GalPerelman https://github.com/GalPerelman https://github.com/GalPerelman Thanks for your reply. I modified PV.adapt(delta) to for t in range(T_len): PV[t].adapt(delta[t]) But the model is still infeasible. Indeed, PV is a LDR with length of 24, and delta is a RVAR with length of 24, denoting the fluctuation of Photovoltaic. The relationship of PV and delta is PV = PV0+delta where PV0 is a given parameter. — Reply to this email directly, view it on GitHub <#46 (comment) <#46 (comment)>>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWDYEQE2N5TKCKBW5MI3XEX5B3ANCNFSM6AAAAAAXX5A7JU. You are receiving this because you are subscribed to this thread. — Reply to this email directly, view it on GitHub <#46 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWD7LQOFJQT4RZJ3YEILXEYAUZANCNFSM6AAAAAAXX5A7JU. You are receiving this because you were mentioned.

chenzjames commented 1 year ago

There is not a direct way in the current version of RSOME. If the uncertain parameter has discrete realizations, then you can treat the probabilities of realizations as the underlying source of uncertainty. The problem might still be a robust optimization problem that could potentially easy to build up in RSOME.

From: Yuxuan @.> Sent: Saturday, May 6, 2023 9:10 PM To: XiongPengNUS/rsome @.> Cc: Zhi Chen @.>; Mention @.> Subject: Re: [XiongPengNUS/rsome] Numerical problems with equality constraints with rvar (Issue #46)

Dr. @chenzjames https://github.com/chenzjames

Thanks a lot, and I wonder how to assign vtype='B' to rvar.

Is it because that RSOME use LDR “approximation”, while CCG yields the exact solution? … On 6 May 2023, at 16:03, Yuxuan @.> wrote: Hi, Dr. @chenzjames https://github.com/chenzjames https://github.com/chenzjames I tried to convert the equality constraints into inequality constraints before, the model is feasible, but the optimization results are far from the results I calculated with the CCG algorithm. How can I avoid the gap? In general, we shall avoid “equality” in a robust constraint, which could be problematic. Most cases, a weak inequality (either >= or <=) will do. … x-msg://6/# On 6 May 2023, at 15:32, Yuxuan @.> wrote: Hi, @GalPerelman https://github.com/GalPerelman https://github.com/GalPerelman https://github.com/GalPerelman Thanks for your reply. I modified PV.adapt(delta) to for t in range(T_len): PV[t].adapt(delta[t]) But the model is still infeasible. Indeed, PV is a LDR with length of 24, and delta is a RVAR with length of 24, denoting the fluctuation of Photovoltaic. The relationship of PV and delta is PV = PV0+delta where PV0 is a given parameter. — Reply to this email directly, view it on GitHub <#46 https://github.com/XiongPengNUS/rsome/issues/46 (comment) <#46 (comment) https://github.com/XiongPengNUS/rsome/issues/46#issuecomment-1537078692 >>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWDYEQE2N5TKCKBW5MI3XEX5B3ANCNFSM6AAAAAAXX5A7JU. You are receiving this because you are subscribed to this thread. — Reply to this email directly, view it on GitHub <#46 (comment) https://github.com/XiongPengNUS/rsome/issues/46#issuecomment-1537084378 >, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWD7LQOFJQT4RZJ3YEILXEYAUZANCNFSM6AAAAAAXX5A7JU. You are receiving this because you were mentioned.

— Reply to this email directly, view it on GitHub https://github.com/XiongPengNUS/rsome/issues/46#issuecomment-1537139322 , or unsubscribe https://github.com/notifications/unsubscribe-auth/AEFJWD33GQVSBMKTKSNVETTXEZESNANCNFSM6AAAAAAXX5A7JU . You are receiving this because you were mentioned. https://github.com/notifications/beacon/AEFJWD3IV3AIFFAA64O7O53XEZESNA5CNFSM6AAAAAAXX5A7JWWGG33NNVSW45C7OR4XAZNMJFZXG5LFINXW23LFNZ2KUY3PNVWWK3TUL5UWJTS3T3RHU.gif Message ID: @. @.> >

sometimesstudy commented 1 year ago

many thanks!