I would like to get the value of the variables at each iteration of the optimization process when using Python GEKKO.
For example, the following code from the documentation (https://gekko.readthedocs.io/en/latest/quick_start.html#example) solves problem HS71:
from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
x = m.Array(m.Var, 4, value=1, lb=1, ub=5)
x1, x2, x3, x4 = x # rename variables
x2.value = 5; x3.value = 5 # change guess
m.Equation(np.prod(x) >= 25) # prod>=25
m.Equation(m.sum([xi**2 for xi in x]) == 40) # sum=40
m.Minimize(x1*x4*(x1 + x2 + x3) + x3) # objective
m.solve()
print(x, m.options.OBJFCNVAL)
The output displayed in the console is the following:
----------------------------------------------------------------
APMonitor, Version 1.0.0
APMonitor Optimization Suite
----------------------------------------------------------------
--------- APM Model Size ------------
Each time step contains
Objects : 1
Constants : 0
Variables : 10
Intermediates: 0
Connections : 5
Equations : 7
Residuals : 7
Number of state variables: 10
Number of total equations: - 7
Number of slack variables: - 1
---------------------------------------
Degrees of freedom : 2
**********************************************
Steady State Optimization with Interior Point Solver
**********************************************
Info: Exact Hessian
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
This is Ipopt version 3.10.2, running with linear solver mumps.
Number of nonzeros in equality constraint Jacobian...: 19
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 10
Total number of variables............................: 10
variables with only lower bounds: 1
variables with lower and upper bounds: 4
variables with only upper bounds: 0
Total number of equality constraints.................: 7
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.6109693e+001 4.00e+001 5.24e-001 0.0 0.00e+000 - 0.00e+000 0.00e+000 0
1 1.6896037e+001 7.28e-001 7.82e-001 -0.9 4.00e+001 - 1.00e+000 1.00e+000h 1
2 1.7121828e+001 1.66e-001 6.51e-001 -1.1 2.15e+000 - 8.23e-001 1.00e+000h 1
3 1.6954231e+001 1.60e-001 6.46e-002 -2.1 6.34e-001 - 1.00e+000 1.00e+000h 1
4 1.7008327e+001 1.68e-002 1.22e-002 -2.9 3.39e-001 - 9.94e-001 1.00e+000h 1
5 1.7013921e+001 3.15e-004 1.47e-004 -4.6 3.56e-002 - 1.00e+000 1.00e+000h 1
6 1.7014017e+001 2.51e-007 4.50e-007 -10.5 5.56e-004 - 9.99e-001 1.00e+000h 1
Number of Iterations....: 6
(scaled) (unscaled)
Objective...............: 1.7014017181012623e+001 1.7014017181012623e+001
Dual infeasibility......: 4.5045555432827566e-007 4.5045555432827566e-007
Constraint violation....: 2.5086595754315365e-007 2.5086595754315365e-007
Complementarity.........: 4.8294565936715572e-008 4.8294565936715572e-008
Overall NLP error.......: 4.5045555432827566e-007 4.5045555432827566e-007
Number of objective function evaluations = 7
Number of objective gradient evaluations = 7
Number of equality constraint evaluations = 7
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 7
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 6
Total CPU secs in IPOPT (w/o function evaluations) = 0.014
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
The solution was found.
The final value of the objective function is 17.014017181012623
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 0.021100000000000008 sec
Objective : 17.014017181012623
Successful solution
---------------------------------------------------
[[1.0000000339] [4.7429996271] [3.8211500196] [1.3794082228]] 17.014017181
This includes the value of the objective function at each iteration, but I would also like to get the value of the variables (x) at each iteration.
The solvers don't have call-backs to Python. One way to get the values at each iteration is to call the solver multiple times with a max_iter (maximum iteration) limit and report the values until the solver is successful. Use debug=0 to not raise an exception with an unsuccessful solution.
from gekko import GEKKO
import numpy as np
max_iter = 0
while True:
m = GEKKO(remote=False)
x = m.Array(m.Var, 4, value=1, lb=1, ub=5)
x1, x2, x3, x4 = x # rename variables
x2.value = 5; x3.value = 5 # change guess
m.Equation(np.prod(x) >= 25) # prod>=25
m.Equation(m.sum([xi**2 for xi in x]) == 40) # sum=40
m.Minimize(x1*x4*(x1 + x2 + x3) + x3) # objective
m.options.MAX_ITER = max_iter
m.solve(disp=False,debug=0)
print(max_iter,x, m.options.OBJFCNVAL)
if m.options.APPSTATUS==1:
# break loop with successful solution
break
else:
# increment maximum iterations
max_iter += 1
Here is the iteration summary with iteration number, variable values x, and objective function value. The objective function may get worse (higher) because the equations are not yet satisfied, and the solution may be infeasible.
0 [[1.00999999] [4.960000049] [4.960000049] [1.00999999]] 16.109692918
1 [[1.1256700056] [4.4664784656] [4.2706611878] [1.1371887857]] 16.896036995
2 [[1.1003048073] [4.6596383632] [3.9630056459] [1.2300025937]] 17.121828401
3 [[1.0211678027] [4.7139849665] [3.8710712025] [1.3337143994]] 16.954231463
4 [[1.0027660895] [4.740529129] [3.8261434464] [1.3737295983]] 17.008327097
5 [[1.0000616664] [4.7429670916] [3.8212258627] [1.3792900699]] 17.013920783
6 [[1.0000000339] [4.7429996271] [3.8211500196] [1.3794082228]] 17.014017181
If you have a large model that you don't want to redefine multiple times, then just put the m.solve() function in a loop with m.options.COLDSTART=1.
from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
x = m.Array(m.Var, 4, value=1, lb=1, ub=5)
x1, x2, x3, x4 = x # rename variables
x2.value = 5; x3.value = 5 # change guess
m.Equation(np.prod(x) >= 25) # prod>=25
m.Equation(m.sum([xi**2 for xi in x]) == 40) # sum=40
m.Minimize(x1*x4*(x1 + x2 + x3) + x3) # objective
max_iter = 0
while True:
m.options.MAX_ITER = max_iter
m.options.COLDSTART=1
m.solve(disp=False,debug=0)
print(max_iter,x, m.options.OBJFCNVAL)
if m.options.APPSTATUS==1:
# break loop with successful solution
break
else:
# increment maximum iterations
max_iter += 1