weixin_43026072
Mr.意昂
2019-05-08 12:16

在不同电脑上使用julia JUMP求解优化问题结果不同?

  • 开发语言

我用julia JUMP求解微电网中的机组组合问题(本质上是一个混合整数线性规划问题),调用CPLEX解法器,优化出的结果(储能功率)有时候会出现无意义的锯齿波,分析目标函数和约束条件,这样的结果是不应该出现的。并且我用两台不同电脑使用julia运行完全相同的代码,优化出的结果也有区别,如下图。有没有人可以帮忙分析一下为什么会出现这样的问题?小弟感激不尽!!!
不同电脑上运行的优化结果

两台电脑上的julia版本均为0.6.4,貌似只有几个package版本不同,如下图:
各台电脑上的package版本

以下为优化问题的代码

## Formulation
# m = Model()
m = Model(solver = CplexSolver())
# Variables
@variables m begin
    P_De[iT = 1:nT],  (start = 0, lowerbound = 0, upperbound =De.Pmax )
    U_De[iT = 1:nT],  (Bin,start = 0)
    M_De[iT = 1:nT],  (Bin,start = 0)

    SOC[iT = 1:nT],  (start = Batt.SOC0, lowerbound = Batt.SOCmin, upperbound = Batt.SOCmax)
    P_ch[iT = 1:nT],  (start = 0, lowerbound = 0, upperbound = Batt.Pmax)
    P_dh[iT = 1:nT],  (start = 0, lowerbound = 0, upperbound = Batt.Pmax)
    U_ch[iT = 1:nT],  (Bin,start = 0)
    U_dh[iT = 1:nT],  (Bin,start = 0)
    M_ch[iT = 1:nT],  (Bin,start = 0)
    M_dh[iT = 1:nT],  (Bin,start = 0)

    P_Grid_b[iT = 1:nT],  (start = 0, lowerbound = 0, upperbound = Pgridmax)
    P_Grid_s[iT = 1:nT],  (start = 0, lowerbound = 0, upperbound = Pgridmax)
    U_Grid_b[iT = 1:nT],  (Bin,start = 0)
    U_Grid_s[iT = 1:nT],  (Bin,start = 0)
end
#Objective Function
@objective(m, Min,
        sum(
            + De.cost_flue * dT/3600 * P_De[iT] + De.cost_start * M_De[iT]
            + Batt.cost_P * dT/3600 * (P_ch[iT] + P_ch[iT]) + Batt.cost_cycle * (M_ch[iT] + M_dh[iT])
            + price_buy[iT] * dT/3600 * P_Grid_b[iT] - price_sell[iT] * dT/3600 * P_Grid_s[iT]
            for iT = 1:nT
        )
)
## Equality Constraints
# 功率平衡约束
@constraint(m, Ceq1[iT = 1:nT],  P_Grid_b[iT] - P_Grid_s[iT]
                             + P_dh[iT] - P_ch[iT]
                             + Ppv[iT]
                             + P_De[iT]
                             - Load[iT] == 0)
# SOC定义约束
@constraint(m, Ceq21[iT = 1], Batt.E * (SOC[1] - Batt.SOC0) == dT/3600 * (P_ch[1] - P_dh[1]))
@constraint(m, Ceq22[iT = 2:nT], Batt.E * (SOC[iT] - SOC[iT-1])  == dT/3600 * (P_ch[iT] - P_dh[iT]))

## Inequality Constraints
# 运行备用约束
@constraint(m, C1[iT = 1:nT], Pgridmax - P_Grid_b[iT] + P_Grid_s[iT] + U_De[iT] * De.Pmax - P_De[iT]
                              >= R * Load[iT])
# 柴油机
@constraint(m, Cd11[iT = 1:nT], P_De[iT] >= U_De[iT] * De.Pmin)
@constraint(m, Cd12[iT = 1:nT], P_De[iT] <= U_De[iT] * De.Pmax)

@constraint(m, Cd21[iT = 2:nT], P_De[iT] - P_De[iT-1] >= dT/3600 * De.dPdown)
@constraint(m, Cd220[iT = 1], P_De[1] <= dT/3600 * De.dPup)
@constraint(m, Cd22[iT = 2:nT], P_De[iT] - P_De[iT-1] <= dT/3600 * De.dPup)

@constraint(m, Cd30[iT = 1], M_De[1] >= U_De[1])
@constraint(m, Cd3[iT = 2:nT], M_De[iT] >= U_De[iT] - U_De[iT-1])

# 储能
@constraint(m, Cb11[iT = 1:nT], P_ch[iT] >= U_ch[iT] * Batt.Pmin)
@constraint(m, Cb12[iT = 1:nT], P_ch[iT] <= U_ch[iT] * Batt.Pmax)
@constraint(m, Cb13[iT = 1:nT], P_dh[iT] >= U_dh[iT] * Batt.Pmin)
@constraint(m, Cb14[iT = 1:nT], P_dh[iT] <= U_dh[iT] * Batt.Pmax)

@constraint(m, Cb210[iT = 1], M_ch[1] >= U_ch[1])
@constraint(m, Cb21[iT = 2:nT], M_ch[iT] >= U_ch[iT] - U_ch[iT-1])
@constraint(m, Cb220[iT = 1], M_dh[1] >= U_dh[1])
@constraint(m, Cb22[iT = 2:nT], M_dh[iT] >= U_dh[iT] - U_dh[iT-1])

@constraint(m, Cb3[iT = 1:nT], U_ch[iT] + U_dh[iT] <= 1)

@constraint(m, Cb41[iT = 2:nT], P_ch[iT] - P_ch[iT-1] >= dT/3600 * Batt.dPdown)
@constraint(m, Cb420[iT = 1], P_ch[1] <= dT/3600 * Batt.dPup)
@constraint(m, Cb42[iT = 2:nT], P_ch[iT] - P_ch[iT-1] <= dT/3600 * Batt.dPup)
@constraint(m, Cb43[iT = 2:nT], P_dh[iT] - P_dh[iT-1] >= dT/3600 * Batt.dPdown)
@constraint(m, Cb440[iT = 1], P_dh[1] <= dT/3600 * Batt.dPup)
@constraint(m, Cb44[iT = 2:nT], P_dh[iT] - P_dh[iT-1] <= dT/3600 * Batt.dPup)

@constraint(m, Cb51[iT = 1:nT], SOC[iT] >= Batt.SOCmin + M_dh[iT] * (Batt.SOCmax - Batt.SOCmin - 0.1))
@constraint(m, Cb52[iT = 1:nT], SOC[iT] <= Batt.SOCmax - M_ch[iT] * (Batt.SOCmax - Batt.SOCmin - 0.1))
# 主网
@constraint(m, Cg11[iT = 1:nT], P_Grid_b[iT] <= U_Grid_b[iT] * Pgridmax)
@constraint(m, Cg12[iT = 1:nT], P_Grid_s[iT] <= U_Grid_s[iT] * Pgridmax)
@constraint(m, Cg13[iT = 1:nT], P_Grid_b[iT] >= U_Grid_b[iT] * Pgridmin)
@constraint(m, Cg14[iT = 1:nT], P_Grid_s[iT] >= U_Grid_s[iT] * Pgridmin)

@constraint(m, Cg2[iT = 1:nT], U_Grid_b[iT] + U_Grid_s[iT] <= 1)


status = solve(m)
  • 点赞
  • 回答
  • 收藏
  • 复制链接分享

0条回答