Nicola Puca
Reputation: 21
Solve a constrained non linear system of equation with optimization - gekko python
i'm very new to Python and I've just discovered Gekko. I'm trying to solve a non linear system of 28 equations using an optimization technique. I would like to treat each of the equations as equality constraints of a problem with a constant objective function. Plus, there are also 17 inequality constraints. Is it a good strategy or not? I've tried to initialize the problem with a random vector x0 of 28 components. But each of the solver (APOPT, IPOPT and BPOPT) gave me a failure.
All the numbers you see in the comments of x0 are the solutions of the system to which I have to tend.
Any advice to proceed?
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()
Upvotes: 2
Views: 361
Answers (1)
John Hedengren
Reputation: 14386
Three suggestions:
- Improve the chance of a successful solution by replacing
abs()withm.abs3()as a version of the absolute value that gives the solver continuous first and second derivatives. - Avoid potential divide-by-zero by reformulating equations:
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)
as equations without the division:
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)
- For infeasible solutions, try commenting out (removing) the constraints until it solves successfully. I removed equations until it solved and then started adding them back again until it failed. The only equation that produces an infeasible solution is:
m.Equation(-x24 + x15 - m.abs3(x25) == 0)
- Add an objective function to guide the values if there are fewer than 28 active constraints (equality constraints + inequality constraints at the bounds). The objective is currently
m.Obj(1000)as a constant.
Here is a version that solves successfully without the one constraint:
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()
Upvotes: 0
Related Questions
- Multiple equations minimized in regression IMODE of gekko
- Solving linear programming by python
- Optimal Control Problem with nonlinear equations
- Solving a non-linear system on Python
- '@error: Solution Not Found' for estimating system of equations
- Formulate a nonlinear optimization problem in Python (gekko) and with m.Equations
- Gekko Solver does not satisfy Constraints gives wrong solution
- Python Gekko - Constrained are getting violated although i have an optimal solution
- Gekko optimization
- Trying to solve this non linear optimization using GEKKO, getting this error