# Optimal control example taken from:
# https://github.com/Pyomo/pyomo/blob/master/examples/dae/Optimal_Control.py
# https://github.com/Pyomo/pyomo/blob/master/examples/dae/run_Optimal_Control.py
#
#
# min X2(tf)
# s.t. X1_dot = u
# X2_dot = X1^2 + u^2
# tf = 1
# X1(0) = 1, X2(0) = 0
#
# set-up python version to use
reticulate::use_python("/opt/anaconda3/bin/python")
library(reticulate)
# import pyomo
pyo = import("pyomo.environ", convert = FALSE)
pyde = import("pyomo.dae", convert = FALSE)
bi = import_builtins()
# had signal handling issues when using solver. Turned
# off signal handling based on this discussion thread
# https://github.com/PyUtilib/pyutilib/issues/31
pyulib = import("pyutilib.subprocess.GlobalData")
pyulib$DEFINE_SIGNAL_HANDLERS_DEFAULT = FALSE
# printing model to a file, reading it back into R to print in R markdown
mdl_print = function(m) {
printstr = "
with open('tmp_model.txt', 'w') as f:
r.mdl.pprint(ostream=f)
f.close()
"
printstr = gsub("mdl", m, printstr)
py_run_string(printstr)
printstr_R = paste(readLines("tmp_model.txt"), collapse = "\n")
return(printstr_R)
}
# get operators
source("operators.R")
# load other libraries
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
# create model object
M = pyo$ConcreteModel()
# set time variable
M$t = pyde$ContinuousSet(bounds = tuple(0, 1))
# set state variables x1, x2
M$x1 = pyo$Var(M$t, bounds = tuple(0, 1))
M$x2 = pyo$Var(M$t, bounds = tuple(0, 1))
# set control variable u
M$u = pyo$Var(M$t, initialize = 0)
# set variables denoting derivatives of state variables x1dot, x2dot
M$x1dot = pyde$DerivativeVar(M$x1, wrt = M$t)
M$x2dot = pyde$DerivativeVar(M$x2, wrt = M$t)
# objective: minimize x2 at tf = 1
M$obj = pyo$Objective(expr = M$x2[1], sense = pyo$minimize)
# define constraint x1dot = u
x1dot_rule = function(m, i) {
if (i == 0) {
return(py_get_attr(pyo$Constraint, "Skip"))
} else {
expr = m$x1dot[i] == m$u[i]
}
return(expr)
}
M$x1dot_cons = pyo$Constraint(M$t, rule = x1dot_rule)
# define constraint x2dot = x1^2 + u^2
x2dot_rule = function(m, i) {
if (i == 0) {
return(py_get_attr(pyo$Constraint, "Skip"))
} else {
expr = m$x2dot[i] == m$x1[i]^2 + m$u[i]^2
}
return(expr)
}
M$x2dot_cons = pyo$Constraint(M$t, rule = x2dot_rule)
# set initial conditions: x1(0) = 1, x2(0) = 0
M$init_x1 = pyo$Constraint(expr = M$x1[0] == 1.0)
M$init_x2 = pyo$Constraint(expr = M$x2[0] == 0.0)
# apply collocation
discretizer = pyo$TransformationFactory('dae.collocation')
discretizer$apply_to(M,nfe=20L,ncp=3L,scheme='LAGRANGE-RADAU')
## None
discretizer$reduce_collocation_points(M,var=M$u,ncp=1L,contset=M$t)
## unknown
solver = pyo$SolverFactory("ipopt")
results = solver$solve(M, logfile = "tmp.log")
cat(paste(readLines("tmp.log"), collapse = "\n"))
## Solver command line: ['/opt/anaconda3/bin/ipopt', '/var/folders/h0/0h0h0yrs6f76syrmv1_q13mw0000gn/T/tmpskotnynz.pyomo.nl', '-AMPL']
##
## Ipopt 3.12.12:
##
## ******************************************************************************
## 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.12, running with linear solver mumps.
## NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
##
## Number of nonzeros in equality constraint Jacobian...: 982
## Number of nonzeros in inequality constraint Jacobian.: 0
## Number of nonzeros in Lagrangian Hessian.............: 120
##
## Total number of variables............................: 302
## variables with only lower bounds: 0
## variables with lower and upper bounds: 122
## variables with only upper bounds: 0
## Total number of equality constraints.................: 282
## 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 9.9999900e-03 9.90e-01 1.39e-04 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
## 1 5.8905427e-01 1.09e+01 1.46e+02 -1.7 1.92e+01 - 5.50e-04 3.21e-02F 1
## 2 4.2017983e-01 2.19e+01 7.14e+03 -1.7 4.65e+00 - 1.40e-04 9.90e-01f 1
## 3 6.8404869e-01 7.34e+00 7.78e+01 -1.7 7.47e+00 - 6.16e-01 9.90e-01h 1
## 4 8.0231348e-01 1.36e+00 3.62e+03 -1.7 1.14e+00 - 8.13e-01 9.90e-01h 1
## 5 8.2251267e-01 1.66e-02 3.17e-02 -1.7 9.48e-01 - 1.00e+00 1.00e+00h 1
## 6 8.1510746e-01 1.18e-02 5.27e+04 -3.8 3.10e-01 - 9.48e-01 1.00e+00h 1
## 7 7.1223952e-01 4.21e-01 1.01e+04 -3.8 1.80e+00 - 8.08e-01 1.00e+00h 1
## 8 7.6169962e-01 8.70e-05 1.37e+03 -3.8 4.12e-01 - 8.65e-01 1.00e+00h 1
## 9 7.6177628e-01 1.58e-04 1.52e-05 -3.8 1.87e-02 - 1.00e+00 1.00e+00h 1
## iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
## 10 7.6165391e-01 6.80e-04 1.27e-05 -5.7 3.92e-02 - 1.00e+00 1.00e+00h 1
## 11 7.6171721e-01 3.58e-07 1.30e-08 -5.7 5.98e-04 - 1.00e+00 1.00e+00h 1
## 12 7.6171723e-01 1.18e-07 2.23e-09 -8.6 5.00e-04 - 1.00e+00 1.00e+00h 1
## 13 7.6171725e-01 4.25e-14 5.15e-16 -8.6 1.00e-07 - 1.00e+00 1.00e+00h 1
##
## Number of Iterations....: 13
##
## (scaled) (unscaled)
## Objective...............: 7.6171724526877993e-01 7.6171724526877993e-01
## Dual infeasibility......: 5.1485543769626019e-16 5.1485543769626019e-16
## Constraint violation....: 2.8227360768296603e-14 4.2521541843143495e-14
## Complementarity.........: 2.5059036763221995e-09 2.5059036763221995e-09
## Overall NLP error.......: 2.5059036763221995e-09 2.5059036763221995e-09
##
##
## Number of objective function evaluations = 15
## Number of objective gradient evaluations = 14
## Number of equality constraint evaluations = 15
## Number of inequality constraint evaluations = 0
## Number of equality constraint Jacobian evaluations = 14
## Number of inequality constraint Jacobian evaluations = 0
## Number of Lagrangian Hessian evaluations = 13
## Total CPU secs in IPOPT (w/o function evaluations) = 0.067
## Total CPU secs in NLP function evaluations = 0.001
##
## EXIT: Optimal Solution Found.
##
# Get the results
#
# NOTE:
# the time points are in the list M$t but it needs to accessed using bi$list(M$t)
# state and control variables are index by elements of M$t.
# However accessing say x1 at time point i using M$x[i]() doesn't work.
# What seems to work is using M$x[tuple(i)]()
#
res_L = list()
k = 1
for(i in bi$list(M$t)) {
res_L[[k]] = list(t = i, x1 = py_to_r(M$x1[tuple(i)]()),
x2 = py_to_r(M$x2[tuple(i)]()),
u = py_to_r(M$u[tuple(i)]()))
k = k + 1
}
res_df = bind_rows(res_L)
# plot the state and control variables
p = ggplot(res_df) + geom_line(aes(x = t, y = x1, color = "x1"))
p = p + geom_line(aes(x = t, y = x2, color = "x2"))
p = p + geom_line(aes(x = t, y = u, color = "u"))
p = p + labs(x = "time", y = "")
p = p + theme_bw()
p
Session Info
sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.2 dplyr_1.0.0 reticulate_1.16
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.4.6 knitr_1.28 magrittr_1.5 munsell_0.5.0
## [5] tidyselect_1.1.0 colorspace_1.4-1 lattice_0.20-41 R6_2.4.1
## [9] rlang_0.4.6 stringr_1.4.0 highr_0.8 tools_4.0.1
## [13] grid_4.0.1 gtable_0.3.0 xfun_0.14 withr_2.2.0
## [17] htmltools_0.5.0 ellipsis_0.3.1 yaml_2.2.1 digest_0.6.25
## [21] tibble_3.0.1 lifecycle_0.2.0 crayon_1.3.4 Matrix_1.2-18
## [25] farver_2.0.3 purrr_0.3.4 vctrs_0.3.1 glue_1.4.1
## [29] evaluate_0.14 rmarkdown_2.3 labeling_0.3 stringi_1.4.6
## [33] compiler_4.0.1 pillar_1.4.4 scales_1.1.1 generics_0.0.2
## [37] jsonlite_1.6.1 pkgconfig_2.0.3
# pyomo version used is 5.7