How to check, how good the Constraint Equation in my Gekko Model solution are fullfilled?

Question

Hello I was using Python Gekko to solve a minimization problem. Currently I am not satisfied with the solution from my Gekko Model and would like to improve the solution by iteratively selecting different start values. Therefore I would like to check how good the Equations/ Constraints in my Gekko Model are fullfilled after each iteration. Is there any solution variable that measures this value? Is it possible to write a generic Utility function that iterates over all equations that are stored in the ._equation variable that is mentioned in the docs and prints out how far the constraint is violated? For example I have the Equation 'a+b==1' and the solution is 'a=1.001' and 'b=0' then my constraint is violated by '0.001'. To print out the results of my solution for each Gekko-variable I wrote a little generic function that takes the Gekko-Model as input and iterates over all variables that are stored in the '._variables' variable and prints out the name of the variable with the numeric value.

I already tried to iterate over the equations and print them out:

for e in m._equations:
    print(e)

It prints out the equation with the variable names inside the equation. For example:

(((2)*(phi_h))+phi_h2)<6.283185307179586

However, what I want is somewhat like replacing all variable names with their current value.

Engagement

Answer 1

There is a built-in report to give details about the equation residuals. For the following sample application, set remote=False and open the run directory with m.open_folder() to find the report rto_4_eqn.txt:

from gekko import GEKKO

m = GEKKO(remote=False)
x = m.Var(value=0,name='x')
y = m.Var(value=1,name='y')
m.Equations([x + 2*y==1, x**2+y**2==1])
m.options.DIAGLEVEL=2
m.solve(disp=False)
print([x.value[0],y.value[0]])
m.open_folder()

If the solver does not find a solution, there is an infeasibilities.txt file that is generated. You can test this by changing the equation to x + 2*y==10 where there is no intersection (solution) between the line and circle. Set debug=0 to avoid throwing an exception if a solution is not found:

from gekko import GEKKO

m = GEKKO(remote=False)
x = m.Var(value=0,name='x')
y = m.Var(value=1,name='y')
m.Equations([x + 2*y==1, x**2+y**2==1])
m.options.DIAGLEVEL=2
m.solve(disp=True,debug=0)
m.open_folder()

If you want to also have this report when remote=True or don't want to open the run directory each time, use the additional code to automatically display the correct file for the equation report whether the solution is successful or not. The file name changes with IMODE value.

from gekko import GEKKO
import os

m = GEKKO(remote=False)
x = m.Var(value=0,name='x')
y = m.Var(value=1,name='y')
m.Equations([x + 2*y==10, x**2+y**2==1])
m.options.DIAGLEVEL=2
m.solve(disp=False,debug=0)
print([x.value[0],y.value[0]])
#m.open_folder()

if m.options.APPSTATUS==1:
    print('Successful')
    im = m.options.IMODE
    if im==1:
        file = 'ss_4_eqn.txt'
    elif im==2:
        file = 'mpu_4_eqn.txt'
    elif im==3:
        file = 'rto_4_eqn.txt'
    elif im==4:
        file = 'sim_4_eqn.txt'
    elif (im==5) or (im==8):
        file = 'est_4_eqn.txt'
    elif (im==6) or (im==9):
        file = 'ctl_4_eqn.txt'
    elif im==7:
        file = 'sqs_4_eqn.txt'
else:
    print('Not successful')
    file = 'infeasibilities.txt'

if m._remote:
    from gekko.apm import get_file
    # remote solve - get file
    f = get_file(m._server,m._model_name,file)
    f = f.decode().replace('\r','')
    with open(file, 'w') as fl:
        fl.write(str(f))
        print(f)
else:
    # local solution - get file from run folder
    with open(os.path.join(m._path,file), 'r') as fl:
        print(fl.read())

There is additional information about the infeasibilities.txt file at How to retrieve the 'infeasibilities.txt' from the gekko

Sample Report infeasibilities.txt with Infeasible Problem

[0.5, 1.0]
Not successful
************************************************
***** POSSIBLE INFEASBILE EQUATIONS ************
************************************************
____________________________________________________________________________
EQ Number   Lower        Residual     Upper        Infeas.     Name
         1   0.0000E+00  -7.5000E+00   0.0000E+00   7.5000E+00  ss.Eqn(1): 0 = (x+((2)*(y)))-(10)
 Variable   Lower        Value        Upper        $Value      Name
         1  -1.2346E+20   5.0000E-01   1.2346E+20   0.0000E+00  ss.x
         2  -1.2346E+20   1.0000E+00   1.2346E+20   0.0000E+00  ss.y
____________________________________________________________________________
EQ Number   Lower        Residual     Upper        Infeas.     Name
         2   0.0000E+00   2.5000E-01   0.0000E+00  -2.5000E-01  ss.Eqn(2): 0 = (((x)^(2))+((y)^(2)))-(1)
 Variable   Lower        Value        Upper        $Value      Name
         1  -1.2346E+20   5.0000E-01   1.2346E+20   0.0000E+00  ss.x
         2  -1.2346E+20   1.0000E+00   1.2346E+20   0.0000E+00  ss.y
************************************************
****** ACTIVE OBJECTIVE EQUATIONS **************
************************************************
Number           ID  Node    Horizon  Unscaled Res  Scaled Res   Scaling     Name
************************************************
************* ACTIVE EQUATIONS *****************
************************************************
Number           ID  Node    Horizon  Unscaled Res  Scaled Res   Scaling     Name
         1          1    1          1  -7.5000E+00  -7.5000E+00   1.0000E+00  ss.Eqn(1): 0 = (x+((2)*(y)))-(10)
         2          2    1          1   2.5000E-01   2.5000E-01   1.0000E+00  ss.Eqn(2): 0 = (((x)^(2))+((y)^(2)))-(1)
************************************************
************ INACTIVE EQUATIONS ****************
************************************************
Number     Unscaled Res    Scaled Res   Scaling      Name

Sample Report rto_4_eqn.txt with Successful Solution

[-0.60000000466, 0.80000000233]
Successful
************************************************
***** POSSIBLE INFEASBILE EQUATIONS ************
************************************************
************************************************
****** ACTIVE OBJECTIVE EQUATIONS **************
************************************************
Number           ID  Node    Horizon  Unscaled Res  Scaled Res   Scaling     Name
************************************************
************* ACTIVE EQUATIONS *****************
************************************************
Number           ID  Node    Horizon  Unscaled Res  Scaled Res   Scaling     Name
         1          1    1          1   0.0000E+00   0.0000E+00   1.0000E+00  ss.Eqn(1): 0 = (x+((2)*(y)))-(1)
         2          2    1          1   9.3234E-09   9.3234E-09   1.0000E+00  ss.Eqn(2): 0 = (((x)^(2))+((y)^(2)))-(1)
************************************************
************ INACTIVE EQUATIONS ****************
************************************************
Number     Unscaled Res    Scaled Res   Scaling      Name

Pregnant

Thanks for your answer. I had already used the 'infeasabilities.txt' file, but your additional informations are very helpful. I am curious: How are the values listed in the 'infeasibilities.txt' file "Residual" and "Infeas." calculated? Whats the difference between both? How are the values in the file 'rto_4_eqn.txt': "Unscaled Res" and "Scaled Res" calculated? Does the "Res" stand for residual? Is it possible to print the results out like: 'm._equations[0].UnscaledRes' or 'm._equations[0].Residual'? – Engagement

Yes, "Res" is the residual. The infeasibility is the non-zero residual when the solver can't find a solution. Gekko applied internal scaling based on the initial guess of the variable. If the value is default then no scaling is applied. I'm not aware a method to calculate m._equations[0].Residual besides performing an eval() function with the variable values (e.g. x.value[0]). – Pregnant