Non-linear optimisation in Python using GEKKO: maximum weight constraint

Question

I'd like to be able to fix an error in the Python code below relating to "wmax". I'm trying to find a solution of weights (the w[i]) such that all the weights are less than or equal to 2.5.

Here is the code:

import numpy as np
import pandas as pd
from gekko import GEKKO
x1=pd.DataFrame([1,5,2,3,2,5,2,6,
                 8,4,9,9,7,5,2,4]).iloc[:,0]
    
m = GEKKO(remote=False)

M1 = 5
TM1 = 0.05
S1=3
TS1=0.05

n = len(x1)
a0,a1,a2,a3,a4 = m.Array(m.Var,5,lb=-500,ub=500)
w = [m.Intermediate(m.exp(a0+a1*x1[i]
                          +a2*x1[i]**2
                          +a3*x1[i]**3
                          +a4*x1[i]**4))
     for i in range(n)]
w2 = [w[i]**2 for i in range(n)]
wx1 = [w[i]*x1[i] for i in range(n)]
wx1M1 = [w[i]*(x1[i]-M1)**2 for i in range(n)]
wmax=  [(((n/m.sum(w))*w[i])<=2.5) for i in range(n)]
m.Equation(m.abs3(m.sum(wx1)/n-M1)<=TM1)
m.Equation(m.abs3(m.sqrt(m.sum(wx1M1)/(n-1))-S1)<=TS1)
m.Equation(m.sum(w)==n)
m.Equation(m.sum(wmax)==n)
m.Minimize(m.sum(w2))
m.options.SOLVER = 1
m.solve()

This is the error I get:

Traceback (most recent call last):
  File "C:\...\sas.py", line 43, in <module>
    m.solve()
  File "C:\...\gekko\gekko.py", line 2140, in solve
 ---------------------------------------------
 APMonitor, Version 1.0.0
 APMonitor Optimization Suite
 ---------------------------------------------
 @error: Duplicate Names
 error: duplicate variable declarations
 intermediate: 1 
 intermediate: 2 
 STOPPING. . .

Answer 1

The problem with wmax is the inequality <= in the list comprehension.

wmax=  [(((n/m.sum(w))*w[i])<=2.5) for i in range(n)]

Instead, separate the definition of wmax as a variable and include an upper bound on the value. The solver attempts to find a solution with this modification, but fails to converge. Try removing the m.abs3() functions:

m.Equation(m.abs3(m.sum(wx1)/n-M1)<=TM1)
m.Equation(m.abs3(m.sqrt(m.sum(wx1M1)/(n-1))-S1)<=TS1)

with an equivalent reformulation:

c1 = m.Var(lb=-TM1,ub=TM1)
m.Equation(c1==m.sum(wx1)/n-M1)

c2 = m.Var(lb=-TS1,ub=TS1)
m.Equation(c2==m.sqrt(m.sum(wx1M1)/(n-1))-S1)

This produces a successful solution in 37 iterations:

 Iter    Objective  Convergence
   30  1.75486E+01  2.10425E-03
   31  1.75001E+01  1.60016E-04
   32  1.74959E+01  9.88379E-07
   33  1.74951E+01  6.02537E-06
   34  1.74948E+01  4.18756E-06
   35  1.74948E+01  1.67332E-07
   36  1.74948E+01  2.26932E-10
   37  1.74948E+01  2.26932E-10
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.11670000000000003 sec
 Objective      :  17.49481149818343
 Successful solution
 ---------------------------------------------------

Here is the complete code:

import numpy as np
import pandas as pd
from gekko import GEKKO
x1=pd.DataFrame([1,5,2,3,2,5,2,6,
                 8,4,9,9,7,5,2,4]).iloc[:,0]
m = GEKKO(remote=False)
M1 = 5; TM1 = 0.05; S1=3; TS1=0.05
n = len(x1)
a0,a1,a2,a3,a4 = m.Array(m.Var,5,lb=-500,ub=500)
w = [m.Intermediate(m.exp(a0+a1*x1[i]
                          +a2*x1[i]**2
                          +a3*x1[i]**3
                          +a4*x1[i]**4))
                       for i in range(n)]
w2 = [w[i]**2 for i in range(n)]
wx1 = [w[i]*x1[i] for i in range(n)]
wx1M1 = [w[i]*(x1[i]-M1)**2 for i in range(n)]
wmax= m.Array(m.Var,n,ub=2.5)
m.Equations([wmax[i]==(n/m.sum(w))*w[i] for i in range(n)])
c1 = m.Var(lb=-TM1,ub=TM1)
m.Equation(c1==m.sum(wx1)/n-M1)
c2 = m.Var(lb=-TS1,ub=TS1)
m.Equation(c2==m.sqrt(m.sum(wx1M1)/(n-1))-S1)
m.Equation(m.sum(w)==n)
m.Equation(m.sum(wmax)==n)
m.Minimize(m.sum(w2))
m.options.SOLVER = 1
m.solve()