BYU-PRISM / GEKKO

GEKKO Python for Machine Learning and Dynamic Optimization
https://machinelearning.byu.edu
Other
580 stars 103 forks source link

unexpected 'Warning: connection failed' when using fix to set initial condition #80

Closed whoburg closed 3 years ago

whoburg commented 4 years ago

This is possibly related to #59, but possibly a separate issue. Here is a GEKKO implementation of a simple optimal control problem -- moving a double integrator (i.e. a frictionless brick) to the origin in a fixed amount of time.

In the code below, fix is used to set both the initial condition (right of the origin and moving right) and the desired final condition of x=0, v=0.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO()
ntime = 51
m.time = np.linspace(0, 2, ntime)

x = m.Var(value=0)
v = m.Var(value=0)
# set initial condition via fix
m.fix(x, pos=0, val=10)  # pos=1 these 2 lines leads to infeasible problem (??)
m.fix(v, pos=0, val=5)

U = m.MV()
U.STATUS = 1
m.Equation(x.dt() == v)
m.Equation(v.dt() == U)

# fix endpoint
m.fix(x, pos=len(m.time) - 1, val=0)
m.fix(v, pos=len(m.time) - 1, val=0)
m.Obj(U**2)
m.options.IMODE = 6
m.solve()

plt.figure(figsize=(4, 3))
plt.subplot(2, 1, 1)
plt.plot(m.time, x.value, label=r"$x$")
plt.plot(m.time, v.value, label=r"$v$")
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(m.time, U.value, label=r"$U$")
plt.legend()
plt.xlabel("Time")

The code runs and generates the expected plots, but also prints four unexpected Warning: connection failed, unable to locate variable messages.

apm 174.240.5.2_gk_model69 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------

 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :            3
   Intermediates:            0
   Connections  :            8
   Equations    :            3
   Residuals    :            3

 Warning: connection failed
 p(0).n(2).v1 = 10
 unable to locate variable
 p(0).n(2).v1
 10

 Warning: connection failed
 p(0).n(2).v1 = fixed
 unable to locate variable
 p(0).n(2).v1
 fixed

 Warning: connection failed
 p(0).n(2).v2 = 5
 unable to locate variable
 p(0).n(2).v2
 5

 Warning: connection failed
 p(0).n(2).v2 = fixed
 unable to locate variable
 p(0).n(2).v2
 fixed

 Number of state variables:            446
 Number of total equations: -          400
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :             46

 **********************************************
 Dynamic Control 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.12.10, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:      740
Number of nonzeros in inequality constraint Jacobian.:      200
Number of nonzeros in Lagrangian Hessian.............:       50

Total number of variables............................:      446
                     variables with only lower bounds:      100
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      300
Total number of inequality constraints...............:      100
        inequality constraints with only lower bounds:      100
   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  3.4499965e-05 1.00e+01 1.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.1852665e-01 9.95e+00 1.53e+00  -0.0 1.43e+02    -  3.45e-03 5.33e-03f  1
   2  1.8412365e+03 5.56e+00 1.94e+01  -0.3 5.74e+02    -  3.92e-03 4.41e-01h  1
   3  3.0575488e+03 4.22e+00 1.47e+01  -1.9 3.90e+02    -  3.23e-01 2.41e-01h  1
   4  9.0530239e+03 1.55e-06 1.52e+00  -0.3 2.83e+02    -  4.37e-01 1.00e+00h  1
   5  9.0498597e+03 2.67e-08 2.74e-01  -0.6 2.64e+01    -  7.30e-01 1.00e+00f  1
   6  9.0497524e+03 3.79e-09 1.40e-02  -2.6 7.33e+00    -  9.31e-01 1.00e+00f  1
   7  9.0497523e+03 3.02e-11 1.65e-04  -4.0 3.45e-01    -  9.88e-01 1.00e+00h  1
   8  9.0497513e+03 8.17e-14 2.27e-13  -9.7 5.15e-01    -  1.00e+00 1.00e+00f  1
   9  9.0497484e+03 2.49e-14 2.74e-06 -10.6 1.94e+00    -  1.00e+00 7.13e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  9.0497479e+03 3.55e-15 6.10e-06  -6.9 4.52e-01    -  1.00e+00 7.69e-01f  1
  11  9.0497479e+03 5.33e-15 3.15e-07 -11.0 7.38e-02    -  1.00e+00 9.49e-01f  1

Number of Iterations....: 11

                                   (scaled)                 (unscaled)
Objective...............:   9.0497479000674884e+03    9.0497479000674884e+03
Dual infeasibility......:   3.1463071921944386e-07    3.1463071921944386e-07
Constraint violation....:   5.3290705182007514e-15    5.3290705182007514e-15
Complementarity.........:   1.2866925199712653e-07    1.2866925199712653e-07
Overall NLP error.......:   2.0480454903723479e-07    3.1463071921944386e-07

Number of objective function evaluations             = 12
Number of objective gradient evaluations             = 12
Number of equality constraint evaluations            = 12
Number of inequality constraint evaluations          = 12
Number of equality constraint Jacobian evaluations   = 12
Number of inequality constraint Jacobian evaluations = 12
Number of Lagrangian Hessian evaluations             = 11
Total CPU secs in IPOPT (w/o function evaluations)   =      0.018
Total CPU secs in NLP function evaluations           =      0.022

EXIT: Optimal Solution Found.

 The solution was found.

 The final value of the objective function is    9049.74790006749

 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   4.780000000027940E-002 sec
 Objective      :    9049.74790006749
 Successful solution
 ---------------------------------------------------

A related question: what is the recommended way to set an initial condition in GEKKO when IMODE=6? Most examples in the documentation seem to set initial conditions using the value kwarg of GEKKO.Var. As a new user I find this somewhat confusing - value seems to be used both to set an initial guess (to then later be manipulated by the optimizer) and to fix initial conditions. Does the behavior really change from a) fixing, for the 0-index of the variable, to b) setting an initial guess, for all other indices of the variable?

APMonitor commented 4 years ago

You set initial conditions when you define your position and velocity.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO()
ntime = 51
m.time = np.linspace(0, 2, ntime)

# set initial conditions [10,5]
x = m.Var(value=10)
v = m.Var(value=5)

U = m.MV()
U.STATUS = 1
m.Equation(x.dt() == v)
m.Equation(v.dt() == U)

# fix endpoint
m.fix(x, pos=len(m.time)-1, val=0)
m.fix(v, pos=len(m.time)-1, val=0)
m.Obj(U**2)
m.options.IMODE = 6
m.open_folder()
m.solve()

plt.figure(figsize=(4, 3))
plt.subplot(2, 1, 1)
plt.plot(m.time, x.value, label=r"$x$")
plt.plot(m.time, v.value, label=r"$v$")
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(m.time, U.value, label=r"$U$")
plt.legend()
plt.xlabel("Time")
plt.show()

Alternatively, you can set x.value = 10 and v.value = 5 before the solve to define the initial conditions. The initial conditions also give a guess for the remainder of the time horizon so that the solver starts from those values for all time points. You can also feed in a vector of values if you'd like to provide a better starting guess. You are correct that value is used for input to the model and to store the results after the solve is complete. You'll also see that x.value is an object with two additional x.valve.value that stores the values and x.value.change that is set to True when the value is updated. This is so that the prior values are used in applications like MPC where the prior solution is time-shifted to initialize the next solution.

The reason that your m.fix() doesn't work for initial conditions is because of 2 reasons:

  1. The underlying APM is 1-based instead of 0-based so the first point is p(1).n(1).x in the collocation horizon. You can see the .apm text file that Gekko writes by using m.open_folder().
  2. The issue that you mentioned about fixing the state value also fixes the derivative. Here is an alternative way to use soft constraints. There is also a way to use hard constraints (without specifying the derivative) to fix any point in the horizon but I'd recommend just using x.value=10 instead for initial conditions.

There are a couple similar applications (pendulum cart and inverted pendulum) on the Machine Learning and Dynamic Optimization Course to yours.

whoburg commented 4 years ago

Thank you for the excellent explanation. Everything above makes sense, but I'm still noticing things that make me think there's a potential bug with the indexing.

  1. Understood that APM is 1-based, but then why do we fix the endpoint to pos=len(m.time)-1 instead of pos=len(m.time)?

  2. Note that if we change two lines in (your updated version of the) example code to

    # fix endpoint
    m.fix(x, pos=len(m.time), val=0)
    m.fix(v, pos=len(m.time), val=0)

    then the code fails with the following traceback

    Traceback (most recent call last):
    File "block-bug-report.py", line 24, in <module>
    m.solve()
    File "...venv/lib/python3.7/site-packages/gekko/gekko.py", line 1965, in solve
    self._write_csv()
    File "...venv/lib/python3.7/site-packages/gekko/gk_write_files.py", line 225, in _write_csv
    t[i[0]+1] = i[1] #index is +1 because of prepended header
    IndexError: index 52 is out of bounds for axis 0 with size 52
  3. In the GEKKO Optimization of Multiple Linked Phases Example, the connection is achieved by the following code:

    # Connect phases together at endpoints
    for i in range(n-1):
    m.Connection(x[i+1],x[i],1,len(m.time)-1,1,nodes)
    m.Connection(x[i+1],'CALCULATED',pos1=1,node1=1)

    Why a connection from node 1 to node N-1 in a 1-indexed setting? E.g. if there are 5 nodes, we connect node 1 to node 4?

  4. In my original example code that awkwardly uses fix to set the initial condition as well as the final condition, I fixed pos=0 and pos=len(m.time)-1, as if the positions were 0-indexed. This code actually produces the expected solution and plots, but also with the Warning: connection failed and unable to locate variable messages. The fact that the code still runs makes me think there might be a strange bug in the conversion from 0-indexing to 1-indexing.

Thanks for entertaining the questions and possible bug report. Amazing software. Really enjoying playing with it.

loganbeal commented 4 years ago

@whoburg These quirks are from the implementation of orthogonal collocation on finite elements. The documentation says that the last node of one finite element overlaps with the first node of the next finite element. Instead of having node[1] on element[n], there is only a reference to the identical point node[-1] on element[n-1]. You can see some of these details in the output files like .t0 and .dxdt with a local solve and m.open_folder().

@APMonitor the documentation for fix and Connection should be clarified.

whoburg commented 4 years ago

Thanks @loganbeal, that explains a lot. So does the pos kwarg of fix refer to an element as opposed to a node? I think that explains a lot -- much of my confusion above was due to thinking pos referred to node numbering as opposed to element numbering. Duh.

On the documentation, https://gekko.readthedocs.io/en/latest/quick_start.html#connections currently says "If pos1 or pos2 is not None, the associated var must be a GEKKO variable and the position is the (0-indexed) time-discretized index of the variable."

I've spent an embarrassing amount of time trying to understand the comment about 0-indexing there. I also think the terminology "index of the variable" is potentially confusing -- perhaps this should refer to either nodes or elements?

APMonitor commented 3 years ago

The next version of Gekko has m.fix_initial(x), m.fix_final(x), m.free_initial(x), and m.free_final(x) to help with these common operations. This will hopefully resolve the concerns about indexing.

APMonitor commented 3 years ago

Another option is to use m.Connection(x[i+1],x[i],1,'end',1,'end'). The 'end' designation is the last node or the last interval in the horizon.

whoburg commented 3 years ago

Excellent, the new helper methods sound like great additions. Looking forward to trying them out.

I still think it would be worth updating the Connection documentation to clarify whether the pos and node arguments refer to node or element indices.

Thanks again for the excellent software and responsiveness.

APMonitor commented 3 years ago

Thanks for the suggestion. I updated the documentation here: https://apmonitor.com/wiki/index.php/Main/OptionApmNodes and will put a link to this in ReadTheDocs.