In [1]:

    

%matplotlib inline
from pymoca.backends.casadi.api import transfer_model
import numpy as np
import casadi as ca
from pylab import plot, hold, legend

model = transfer_model('../models/', 'Exponential', {})
dae = model.dae_residual_function
X0 = model.states[0].start

# Discretization steps
N = 10

# Discretization time step
dt = 0.1

# Optimization variables and bounds
X = ca.MX.sym('X', N)
lbX = np.full(N, -np.inf)
ubX = np.full(N, np.inf)

U = ca.MX.sym('U', N)
lbU = np.full(N, 0)
ubU = np.full(N, 1.75)

# Collocate DAE using backwards Euler method
g = []
i = 0
res = dae(0, X[0], (X[0] - X0) / dt, U[0], U[i], ca.MX(), ca.MX())
g.append(res)
for i in range(1, N):
    res = dae(i * dt, X[i], (X[i] - X[i - 1]) / dt, ca.MX(), U[i], ca.MX(), ca.MX())
    g.append(res)
    
# Optimization objective
f = (X[-1] - 1.0)**2

# Solve NLP
nlp = {'x': ca.vertcat(X, U), 'f': f, 'g': ca.vertcat(*g)}
ipopt_options = {'tol': 1e-6}
solver = ca.nlpsol('nlp', 'ipopt', nlp, {'ipopt': ipopt_options})
solution = solver(lbx=ca.vertcat(lbX, lbU), ubx=ca.vertcat(ubX, ubU), lbg=0, ubg=0)

# Plot solution
T = np.linspace(0, N * dt, N + 1)
X = solution['x'][0:N]
U = solution['x'][N:2*N]
plot(T, ca.vertcat(X0, X), label="$x$")
plot(T[1:], U, label="$u$")
legend()


    

    




******************************************************************************
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.11.9, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:       29
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        1

Total number of variables............................:       20
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       10
                     variables with only upper bounds:        0
Total number of equality constraints.................:       10
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.0000000e+00 1.00e-02 1.75e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  8.6929583e-01 5.55e-17 1.03e-02  -1.0 1.01e-01    -  9.39e-01 1.00e+00f  1
   2  3.3357878e-01 4.44e-16 1.62e-02  -1.7 1.41e+00    -  7.02e-01 1.00e+00f  1
   3  1.0793001e-01 4.44e-16 2.11e-02  -2.5 6.11e-01    -  6.56e-01 1.00e+00f  1
   4  3.1680767e-02 8.88e-16 1.11e-16  -2.5 4.92e-01    -  1.00e+00 1.00e+00f  1
   5  6.8183109e-03 8.88e-16 2.43e-03  -3.8 2.15e-01    -  8.60e-01 1.00e+00f  1
   6  1.1418229e-03 6.66e-16 5.55e-17  -3.8 1.49e-01    -  1.00e+00 1.00e+00f  1
   7  6.7803447e-05 5.55e-16 1.73e-04  -5.7 7.15e-02    -  9.63e-01 1.00e+00f  1
   8  8.2258868e-07 7.77e-16 3.14e-17  -5.7 2.15e-02    -  1.00e+00 1.00e+00f  1
   9  2.6202735e-10 5.55e-16 6.10e-17  -7.0 2.74e-03    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 9

                                   (scaled)                 (unscaled)
Objective...............:   2.6202735175817607e-10    2.6202735175817607e-10
Dual infeasibility......:   6.1040582310933900e-17    6.1040582310933900e-17
Constraint violation....:   5.5511151231257827e-16    5.5511151231257827e-16
Complementarity.........:   2.8411303641615059e-07    2.8411303641615059e-07
Overall NLP error.......:   2.8411303641615059e-07    2.8411303641615059e-07


Number of objective function evaluations             = 10
Number of objective gradient evaluations             = 10
Number of equality constraint evaluations            = 10
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 10
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 9
Total CPU secs in IPOPT (w/o function evaluations)   =      0.005
Total CPU secs in NLP function evaluations           =      0.002

EXIT: Optimal Solution Found.
               t_proc [s]   t_wall [s]    n_eval
         nlp      0.00739      0.00688         1
       nlp_f      2.5e-05     2.48e-05        10
       nlp_g     0.000253     0.000253        10
  nlp_grad_f      6.4e-05     6.21e-05        11
  nlp_hess_l     0.000489      0.00049         9
   nlp_jac_g      0.00099      0.00099        11






    
Out[1]:





<matplotlib.legend.Legend at 0x7fe76c6fcf28>