AMPLpy电池套利利润最大化LP编程问题求解求助

问题描述 #

我尝试在24小时周期内，通过控制电池充电功率r_i或放电功率q_i，结合批发电价P_elec,i最大化套利利润（η为电池效率，z为荷电状态，用于限制放电功率q）。设置了包含高低电价的测试用例，理论上应该有明确的利润最大化操作，但求解结果中大部分r_i和q_i为0，仅q[24]有值。

AMPLpy代码 #

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # 安装的模块
    license_uuid="default",  # 使用的许可证
)  # 实例化AMPL对象并注册魔法命令
%%ampl_eval
#define parameters
reset ; 
param R_cap;      # 峰值设备额定功率容量（MWh）
param D;            # 电池储能系统持续时长（小时）
param eta;  # 电池往返效率
param n;
#param P_batt;            # 电池每千瓦价格
param z{i in 1..n};      # 时刻i的荷电状态（kWh）
#param Q{i in 1..n};      # 时刻i的放电功率（kW）
param P_elec{i in 1..n}; # 时刻i的电价（$/kWh）
# define variables

#var R;                    
var r{i in 1..n};        # 时刻i的充电量（kWh）
var q{i in 1..n};        # 时刻i的放电量

# define model and constraints
maximize Profit:
sum{i in 1..n}(P_elec[i]*(q[i]-r[i])); # 最大化套利利润(1)
subject to charging_constraint{i in 1..n}: 0<=r[i]<=R_cap*D-z[i];
subject to discharging_constraint{i in 1..n}: 0<=q[i]<=z[i];       
subject to stat_of_charge{i in 1..n}: 0<=z[i]<=R_cap*D;
subject to charge_discharge_balance{i in 2..n}:
z[i]=z[i-1]+eta*r[i-1]-(1/eta)*q[i-1];             
subject to charge_discharge{i in 1..n}: r[i]*q[i] = 0;  
# load parameter data
ampl.param["R_cap"] = 498.0
#ampl.param["P_batt"]=400.0
ampl.param["D"]=4.0;
ampl.param["n"]=24; 
ampl.param["eta"]=(0.85)**(1/2);
#param["z"]=[0.5*ampl.param["R_cap"]*ampl.param["D"],0] + [0]*23
ampl.param["z"]=[0.5*498.0*4.0]*24
ampl.param["P_elec"]=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,100,100,100,100,1,1,1,1];
ampl.solve(solver="highs")

求解输出 #

HiGHS 1.11.0: 
------------ WARNINGS ------------
WARNING.  72 case(s) of "PLApproxDomain_PowConstExp". One of them:
  Argument domain of a 'PowConstExp'
has been reduced to [0.000000, 316.227766] for numerical reasons
(partially controlled by cvt:plapprox:domain.)
WARNING.  72 case(s) of "PLApprox_PowConstExp". One of them:
  An expression of type 'PowConstExp' has been
piecewise-linearly approximated. Set cvt:plapprox:reltol
to control precision (currently 0.010000).
HiGHS 1.11.0: optimal solution; objective 316.227766
162 simplex iterations
1 branching nodes
 
------------ WARNINGS ------------
WARNING.  72 case(s) of "PLApproxDomain_PowConstExp". One of them:
  Argument domain of a 'PowConstExp'
has been reduced to [0.000000, 316.227766] for numerical reasons
(partially controlled by cvt:plapprox:domain.)
WARNING.  72 case(s) of "PLApprox_PowConstExp". One of them:
  An expression of type 'PowConstExp' has been
piecewise-linearly approximated. Set cvt:plapprox:reltol
to control precision (currently 0.010000).

变量解结果 #

solution = ampl.get_solution(zeros=True)
result = ampl.get_value("solve_result")
print(f"Solve result: {result}")

执行后输出：

{'r[1]': 0, 'r[2]': 0, 'r[3]': 0, 'r[4]': 0, 'r[5]': 0, 'r[6]': 0, 'r[7]': 0, 'r[8]': 0, 'r[9]': 0, 'r[10]': 0, 'r[11]': 0, 'r[12]': 0, 'r[13]': 0, 'r[14]': 0, 'r[15]': 0, 'r[16]': 0, 'r[17]': 0, 'r[18]': 0, 'r[19]': 0, 'r[20]': 0, 'r[21]': 0, 'r[22]': 0, 'r[23]': 0, 'r[24]': 0, 'q[1]': 0, 'q[2]': 0, 'q[3]': 0, 'q[4]': 0, 'q[5]': 0, 'q[6]': 0, 'q[7]': 0, 'q[8]': 0, 'q[9]': 0, 'q[10]': 0, 'q[11]': 0, 'q[12]': 0, 'q[13]': 0, 'q[14]': 0, 'q[15]': 0, 'q[16]': 0, 'q[17]': 0, 'q[18]': 0, 'q[19]': 0, 'q[20]': 0, 'q[21]': 0, 'q[22]': 0, 'q[23]': 0, 'q[24]': 316.227766016838}

荷电状态z结果 #

%%ampl_eval
display z;

执行后输出：

z [*] :=
 1 996    4 996    7 996   10 996   13 996   16 996   19 996   22 996
 2 996    5 996    8 996   11 996   14 996   17 996   20 996   23 996
 3 996    6 996    9 996   12 996   15 996   18 996   21 996   24 996
;

问题分析与解决方案 #

核心问题 #

1. 荷电状态z被错误定义为参数：当前代码中z是参数而非变量，导致荷电状态无法随充放电操作动态变化，所有时刻的z固定为初始值996，约束条件失去作用。
2. 充放电互斥约束非线性：r[i]*q[i] = 0是二次约束，HiGHS作为线性规划求解器只能通过近似处理，可能导致求解异常。
3. 初始状态约束缺失：未明确设置初始荷电状态的约束逻辑，且未考虑周期结束时的荷电状态要求（如回到初始值）。

修正步骤 #

1. 将z改为变量：删除param z{i in 1..n};，添加var z{i in 1..n} >=0 <= R_cap*D;，同时移除原stat_of_charge约束（已包含在变量定义中）。
2. 替换非线性充放电互斥约束：改用线性约束r[i] <= M*(1 - b[i])和q[i] <= M*b[i]，其中b[i]是0-1变量，M为足够大的常数（如电池最大容量）。
3. 明确初始荷电状态：添加约束z[1] = 0.5*R_cap*D;，若需要周期闭环，可添加z[24] = z[1];。
4. 验证单位一致性：确认r、q的单位（kWh）与电价（$/kWh）的匹配性，确保目标函数逻辑正确。

修正后的核心代码片段 #

# 重新定义变量
var r{i in 1..n} >=0 <= R_cap*D;        # 时刻i的充电量（kWh）
var q{i in 1..n} >=0 <= R_cap*D;        # 时刻i的放电量
var z{i in 1..n} >=0 <= R_cap*D;        # 时刻i的荷电状态（kWh）
var b{i in 1..n} binary;                # 充放电互斥二进制变量

# 修正约束
subject to charging_constraint{i in 1..n}: r[i] <= R_cap*D - z[i];
subject to discharging_constraint{i in 1..n}: q[i] <= z[i];       
subject to charge_discharge_balance{i in 2..n}:
z[i] = z[i-1] + eta*r[i-1] - (1/eta)*q[i-1];             
subject to charge_exclusion{i in 1..n}: r[i] <= (R_cap*D)*(1 - b[i]);
subject to discharge_exclusion{i in 1..n}: q[i] <= (R_cap*D)*b[i];
subject to initial_soc: z[1] = 0.5*R_cap*D;

内容的提问来源于stack exchange，提问作者Guilherme Larangeira