用壁虎蟒求解一个受约束的非线性方程组

Question

我对Python很陌生，我刚刚发现了Gekko。 我正在尝试使用优化技术求解一个由 28 个方程组成的非线性系统。我想将每个方程视为具有常数目标函数的问题的相等约束。此外，还有17个不平等约束。 这是一个好的策略吗？ 我尝试使用包含 28 个分量的随机向量 x0 初始化问题。 但是每个求解器(APOPT，IPOPT和BPOPT)都给了我一个失败。

您在x0的评论中看到的所有数字都是我必须倾向于的系统的解决方案。

有什么建议吗？

from gekko import GEKKO
import random
m = GEKKO()
x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28 = [
m.Var() for i in range(28)]
x0 = []
for i in range(29):
x0.append(random.random())
x1.value = x0[1]  # 14.1
x2.value = x0[2]  # 1000
x3.value = x0[3]  # 1000
x4.value = x0[4]  # 1000
x5.value = x0[5]  # 8.1
x6.value = x0[6]  # 960
x7.value = x0[7]  # 965
x8.value = x0[8]  # 8.9
x9.value = x0[9]  # 100.1
x10.value = x0[10]  # 43.1
x11.value = x0[11]  # 23.5
x12.value = x0[12]  # 0.008
x13.value = x0[13]  # 0.543
x14.value = x0[14]  # 0.310
x15.value = x0[15]  # 0.157
x16.value = x0[16]  # 0.106
x17.value = x0[17]  # 0.107
x18.value = x0[18]  # 0.637
x19.value = x0[19]  # 0.845
x20.value = x0[20]  # 1.002
x21.value = x0[21]  # 0.894
x22.value = x0[22]  # 0.893
x23.value = x0[23]  # 0.363
x24.value = x0[24]  # 0.155
x25.value = x0[25]  # -0.002
x26.value = x0[26]  # 0.007
x27.value = x0[27]  # 0.013
x28.value = x0[28]  # 0.102
Ua = m.Param(value=1000)
Ja = m.Param(value=1)
IE = m.Param(value=12.1)
CE = m.Param(value=0.25)
CI = m.Param(value=0.07)
CT = m.Param(value=1 - CE - CI)
p1 = m.Param(value=0.06)
p2 = m.Param(value=0.119)
p3 = m.Param(value=0.160)
p4 = m.Param(value=0.254)
phi0 = m.Param(value=0)
T0 = m.Param(value=1.98)
p0 = m.Param(value=0.002)
m.Obj(1000)
m.Equation(-Ua * Ja + (x24 * x4 + (x23 - x24) * x3 + (x22 - x23) * x2 + (x21 - x22) * x1) + (IE * (x12 + x13 + x14 + x15)) + (
CE * x16 * (1 - p1) * x1 + CE * x17 * (1 - p2) * (x2 - x1 + x8) + CE * x18 * (1 - p3) * (x3 - x2 + x9) + CE * x19 * (1 - p4) * (x4 - x3 + x10)) + (
x16 * p1 * (x5 - phi0 + x1 - x5 + IE) + x17 * p2 * (x6 - x1 + x2 - x7 + IE + x8) + x18 * p3 * (x8 - x2 + x3 - x7 + IE + x9)) + (
x19 * p4 * (Ua - x3 + x10) + (x15 + x19 * (1 - p4)) * (x4 - Ua + x11) - x25 * (x4 - Ua)) == 0)
m.Equation(-x16 * (1 - p1) * (x1 - phi0 + T0) - x16 * p1 * (x5 - phi0 + T0) + x16 * p1 * (x5 - phi0 + T0) + (
x16 * (1 - p1) + x12) * x9 + x12 * IE + x16 * (1 - p1) * CE * (x1 - phi0 + T0) == 0)
m.Equation(-x17 * (1 - p2) * (x2 - x1 + x8) - x17 * p2 * (x6 - x1 + x8) + x17 * p2 * (x6 - x1 + x8) + (
x17 * (1 - p2) + x13) * x9 + x13 * IE + x17 * (1 - p2) * CE * (x2 - x1 + x8) == 0)
m.Equation(-x18 * (1 - p3) * (x3 - x2 + x9) - x18 * p3 * (x7 - x2 + x9) + x18 * p3 * (x7 - x2 + x9) + (
x18 * (1 - p3) + x14) * x10 + x14 * IE + x18 * (1 - p3) * CE * (x3 - x2 + x9) == 0)
m.Equation(-x19 * (1 - p4) * (x4 - x3 + x10) - x19 * p4 * (Ua - x3 + x10) + x19 * p4 * (Ua - x3 + x11) + (
x19 * (1 - p4) + x15) * x11 + x15 * IE + x19 * (1 - p4) * CE * (x4 - x3 + x10) == 0)
m.Equation(-Ja + x20 + x25 == 0)
m.Equation(-Ja + x19 + x24 == 0)
m.Equation(-Ja + x18 + x23 == 0)
m.Equation(-Ja + x17 + x22 == 0)
m.Equation( -Ja + x16 + x21 == 0)
m.Equation(-x16 * p1 + x26 == 0)
m.Equation(-x17 * p2 + x27 == 0)
m.Equation(-x18 * p3 + x28 == 0)
m.Equation(-x16 + p0 * x1 ** (3 / 2) == 0)
m.Equation(-x17 + x16 * (1 - p1) + x12 == 0)
m.Equation(-x18 + x17 * (1 - p2) + x13 == 0)
m.Equation(-x19 + x18 * (1 - p3) + x14 == 0)
m.Equation(-x12 + x16 * (1 - p1) * CI * ((x1 - phi0 + T0) / IE) == 0)
m.Equation(-x13 + x17 * (1 - p2) * CI * ((x2 - x1 + x8) / IE) == 0)
m.Equation(-x14 + x18 * (1 - p3) * CI * ((x3 - x2 + x9) / IE) == 0)
m.Equation(-x15 + x19 * (1 - p4) * CI * ((x4 - x3 + x10) / IE) == 0)
m.Equation(-x21 + x22 + x12 - x26 == 0)
m.Equation(-x22 + x23 + x13 - x27 == 0)
m.Equation(-x23 + x24 + x14 - x28 == 0)
m.Equation(-x24 + x15 - abs(x25) == 0)
m.Equation(-x9 + (CT * x17 * (1 - p2) * (x2 - x1 + x8)) / x18 == 0)
m.Equation(-x10 + (CT * x18 * (1 - p3) * (x3 - x2 + x9)) / x19 == 0)
m.Equation(-x11 + (CT * x19 * (1 - p4) * (x4 - x3 + x10)) / (x19 * (1 - p4) + x15) == 0)
m.Equation(-x8 < 0)
m.Equation(-x9 < 0)
m.Equation(-x10 < 0)
m.Equation(-x11 < 0)
m.Equation(x5 - x1 < 0)
m.Equation(x6 - x2 < 0)
m.Equation(x7 - x3 < 0)
m.Equation(-abs(x24) + abs(x25) < 0)
m.Equation(-x1 < 0)
m.Equation(x1 - x2 < 0)
m.Equation(x2 - x3 < 0)
m.Equation(x3 - x4 < 0)
m.Equation(-x4 < -Ua)
m.Equation(-x5 < 0)
m.Equation(x5 - x6 < 0)
m.Equation(x6 - x7 < 0)
m.Equation(x25 < 0)
m.options.MAX_ITER = 10000
m.options.SOLVER = 0
m.options.IMODE = 2
m.options.DIAGLEVEL = 1
# m.options.DBS_WRITE = 0
# m.options.DBS_READ = 0
# m.options.DBS_LEVEL = 0
m.options.RTOL = 1.0e-3
m.open_folder()
m.solve()

Answer 1

三个建议：

1. 通过将abs()替换为m.abs3()作为绝对值的版本，为求解器提供连续的一阶和二阶导数，从而提高成功求解的机会。
2. 通过重新制定方程来避免潜在的除以零：

m.Equation(-x9 + (CT * x17 * (1 - p2) * (x2 - x1 + x8)) / x18 == 0)
m.Equation(-x10 + (CT * x18 * (1 - p3) * (x3 - x2 + x9)) / x19 == 0)
m.Equation(-x11 + (CT * x19 * (1 - p4) * (x4 - x3 + x10)) / (x19 * (1 - p4) + x15) == 0)

作为没有除法的方程：

m.Equation(-x9 * x18 + (CT * x17 * (1 - p2) * (x2 - x1 + x8)) == 0)
m.Equation(-x10 * x19 + (CT * x18 * (1 - p3) * (x3 - x2 + x9)) == 0)
m.Equation(-x11 * (x19 * (1 - p4) + x15) + (CT * x19 * (1 - p4) * (x4 - x3 + x10)) == 0)

3. 对于不可行的解决方案，请尝试注释掉(删除)约束，直到成功解决。我删除了方程，直到它求解，然后又开始将它们加回去，直到它失败。产生不可行解的唯一方程是：

m.Equation(-x24 + x15 - m.abs3(x25) == 0)

4. 添加一个目标函数来指导值，如果活动约束少于 28 个(边界处的相等约束 + 不等式约束)。目标目前m.Obj(1000)为常量。

这是一个在没有一个约束的情况下成功求解的版本：

from gekko import GEKKO
import random
m = GEKKO()
x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28 = [
m.Var() for i in range(28)]
x0 = []
for i in range(29):
x0.append(random.random())
x1.value = x0[1]  # 14.1
x2.value = x0[2]  # 1000
x3.value = x0[3]  # 1000
x4.value = x0[4]  # 1000
x5.value = x0[5]  # 8.1
x6.value = x0[6]  # 960
x7.value = x0[7]  # 965
x8.value = x0[8]  # 8.9
x9.value = x0[9]  # 100.1
x10.value = x0[10]  # 43.1
x11.value = x0[11]  # 23.5
x12.value = x0[12]  # 0.008
x13.value = x0[13]  # 0.543
x14.value = x0[14]  # 0.310
x15.value = x0[15]  # 0.157
x16.value = x0[16]  # 0.106
x17.value = x0[17]  # 0.107
x18.value = x0[18]  # 0.637
x19.value = x0[19]  # 0.845
x20.value = x0[20]  # 1.002
x21.value = x0[21]  # 0.894
x22.value = x0[22]  # 0.893
x23.value = x0[23]  # 0.363
x24.value = x0[24]  # 0.155
x25.value = x0[25]  # -0.002
x26.value = x0[26]  # 0.007
x27.value = x0[27]  # 0.013
x28.value = x0[28]  # 0.102
Ua = m.Param(value=1000)
Ja = m.Param(value=1)
IE = m.Param(value=12.1)
CE = m.Param(value=0.25)
CI = m.Param(value=0.07)
CT = m.Param(value=1 - CE - CI)
p1 = m.Param(value=0.06)
p2 = m.Param(value=0.119)
p3 = m.Param(value=0.160)
p4 = m.Param(value=0.254)
phi0 = m.Param(value=0)
T0 = m.Param(value=1.98)
p0 = m.Param(value=0.002)
m.Obj(1000)
m.Equation(-Ua * Ja + (x24 * x4 + (x23 - x24) * x3 + (x22 - x23) * x2 + (x21 - x22) * x1) + (IE * (x12 + x13 + x14 + x15)) + (
CE * x16 * (1 - p1) * x1 + CE * x17 * (1 - p2) * (x2 - x1 + x8) + CE * x18 * (1 - p3) * (x3 - x2 + x9) + CE * x19 * (1 - p4) * (x4 - x3 + x10)) + (
x16 * p1 * (x5 - phi0 + x1 - x5 + IE) + x17 * p2 * (x6 - x1 + x2 - x7 + IE + x8) + x18 * p3 * (x8 - x2 + x3 - x7 + IE + x9)) + (
x19 * p4 * (Ua - x3 + x10) + (x15 + x19 * (1 - p4)) * (x4 - Ua + x11) - x25 * (x4 - Ua)) == 0)
m.Equation(-x16 * (1 - p1) * (x1 - phi0 + T0) - x16 * p1 * (x5 - phi0 + T0) + x16 * p1 * (x5 - phi0 + T0) + (
x16 * (1 - p1) + x12) * x9 + x12 * IE + x16 * (1 - p1) * CE * (x1 - phi0 + T0) == 0)
m.Equation(-x17 * (1 - p2) * (x2 - x1 + x8) - x17 * p2 * (x6 - x1 + x8) + x17 * p2 * (x6 - x1 + x8) + (
x17 * (1 - p2) + x13) * x9 + x13 * IE + x17 * (1 - p2) * CE * (x2 - x1 + x8) == 0)
m.Equation(-x18 * (1 - p3) * (x3 - x2 + x9) - x18 * p3 * (x7 - x2 + x9) + x18 * p3 * (x7 - x2 + x9) + (
x18 * (1 - p3) + x14) * x10 + x14 * IE + x18 * (1 - p3) * CE * (x3 - x2 + x9) == 0)
m.Equation(-x19 * (1 - p4) * (x4 - x3 + x10) - x19 * p4 * (Ua - x3 + x10) + x19 * p4 * (Ua - x3 + x11) + (
x19 * (1 - p4) + x15) * x11 + x15 * IE + x19 * (1 - p4) * CE * (x4 - x3 + x10) == 0)
m.Equation(-Ja + x20 + x25 == 0)
m.Equation(-Ja + x19 + x24 == 0)
m.Equation(-Ja + x18 + x23 == 0)
m.Equation(-Ja + x17 + x22 == 0)
m.Equation( -Ja + x16 + x21 == 0)
m.Equation(-x16 * p1 + x26 == 0)
m.Equation(-x17 * p2 + x27 == 0)
m.Equation(-x18 * p3 + x28 == 0)
m.Equation(-x16 + p0 * x1** (3 / 2) == 0)
m.Equation(-x17 + x16 * (1 - p1) + x12 == 0)
m.Equation(-x18 + x17 * (1 - p2) + x13 == 0)
m.Equation(-x19 + x18 * (1 - p3) + x14 == 0)
m.Equation(-x12 + x16 * (1 - p1) * CI * ((x1 - phi0 + T0) / IE) == 0)
m.Equation(-x13 + x17 * (1 - p2) * CI * ((x2 - x1 + x8) / IE) == 0)
m.Equation(-x14 + x18 * (1 - p3) * CI * ((x3 - x2 + x9) / IE) == 0)
m.Equation(-x15 + x19 * (1 - p4) * CI * ((x4 - x3 + x10) / IE) == 0)
m.Equation(-x21 + x22 + x12 - x26 == 0)
m.Equation(-x22 + x23 + x13 - x27 == 0)
m.Equation(-x23 + x24 + x14 - x28 == 0)
#m.Equation(-x24 + x15 - m.abs3(x25) == 0)
m.Equation(-x9 * x18 + (CT * x17 * (1 - p2) * (x2 - x1 + x8)) == 0)
m.Equation(-x10 * x19 + (CT * x18 * (1 - p3) * (x3 - x2 + x9)) == 0)
m.Equation(-x11 * (x19 * (1 - p4) + x15) + (CT * x19 * (1 - p4) * (x4 - x3 + x10)) == 0)
m.Equation(-x8 < 0)
m.Equation(-x9 < 0)
m.Equation(-x10 < 0)
m.Equation(-x11 < 0)
m.Equation(x5 - x1 < 0)
m.Equation(x6 - x2 < 0)
m.Equation(x7 - x3 < 0)
m.Equation(-m.abs3(x24) + m.abs3(x25) < 0)
m.Equation(-x1 < 0)
m.Equation(x1 - x2 < 0)
m.Equation(x2 - x3 < 0)
m.Equation(x3 - x4 < 0)
m.Equation(-x4 < -Ua)
m.Equation(-x5 < 0)
m.Equation(x5 - x6 < 0)
m.Equation(x6 - x7 < 0)
m.Equation(x25 < 0)
m.options.MAX_ITER = 10000
m.options.SOLVER = 1
m.options.IMODE = 3
# m.options.DIAGLEVEL = 1
# m.options.DBS_WRITE = 0
# m.options.DBS_READ = 0
# m.options.DBS_LEVEL = 0
# m.options.RTOL = 1.0e-3
# m.open_folder()
m.solve()