Van der Pol Oscillator

State variables: $y_1(t)$ $y_2(t)$

Control variables: $u(t)$.

Dynamics: $\frac{d y_1}{dt} = y_2$ $\frac{d y_2}{dt} = (1 - {y_1}^2)y_2 - y_1 + u$

Outputs: No outputs for this problem Objective: Minimize:

State Constraints: For this probelem there are no state constraints

Input Constraints: For this problem there are no input constraints

Path Constrains: For this problem there are no path constraints

Integral Constraints: For this problem there are no integral constraints

Start code Include the necessary packages. JuMP is required to setup various configurations while Ipopt is the solver to be used. Technically all other nonlinear solvers available throught JuMP can be used but those have not yet been tested.

# include("../src/DirectOptimalControl.jl")
# import .DirectOptimalControl as DOC

import DirectOptimalControl as DOC

using JuMP
import Ipopt

Let us set first create an optimal control problem. The structure which stores all the data related to the optimal control problem is called OCP (Optimal control problem).

OC = DOC.OCP()

Now we will set various parameters for the solver

OC.tol : Sets up tolerence for the solver. In the intial run it is advisable to keep the tolerance high.
OC.meshitermax : This is the maximum number of iterations that the solver takes
OC.objective_sense: You can set two options here "Max" or "Min" depending on weather the objective is to be minimized or maximized

OC.tol = 1e-7
OC.mesh_iter_max = 20
OC.objective_sense = "Min"

Now it is time to select an optimizer. Let us select the Ipopt optimizer. The OC struct contains an JuMP model in its field OC.model. We can assign any non-linear optimizer that JuMP supports. For this reason we needed to import JuMP. It is advisable to initially keep the solver tolerance higher so that the optimizer converges. You can set all the solver specific options using the JuMP interface to aid the convergence of the solver.

set_optimizer(OC.model, Ipopt.Optimizer)
set_attribute(OC.model, "print_level", 0)
# set_attribute(OC.model, "max_iter", 500)
# set_attribute(OC.model, "tol", 1e-4)

Each OCP must contain atleast one phase. The synatax for adding the phase to an OCP is given by

ph = DOC.PH(OC)

Define the models and cost functions

Now let us define the parameters and functions which make up the model

y10 = 1                      # Initial height
y20 = 0                      # Initial velocity
p1 = -0.4

x0 = [y10, y20]           # Initial state

ph.nk = 0                               # Number of auxillary phase parameters to be optimized
ph.k = @variable(OC.model, [1:ph.nk])   # Assigne them to field k in struct `ph`
OC.nkg = 0                              # Number of auxillary global parameters to be optimized
OC.kg = @variable(OC.model, [1:OC.nkg]) # Assign them to field k in struct `OC`

Auxillary parameters

Now we create a named tuple of various parameters which will be necessary while defining the problem Note that in addition to the constants defined above we can also pass two additional parameters PH.kp and OC.kg. These are additional JuMP variables which can be optimized if required. They will not be used in this example.

p = (p1 = p1, x0 = x0, kp = ph.k, kg = OC.kg )

ns = 2
nu = 1
n = 20

System dynamics

Note that the dyn function which defines the dynamics must be in a particular format. It must five inputs: x : The state of system at time t u : The input of system at time t t : The time t p : p is a named tuple of auzillary parameters required to define the fucntion

function dyn(x, u, t, p)
    y1 = x[1]; y2 = x[2];
    u = u[1]

    y1n = y2
    y2n = (1-y1^2)*y2 - y1 + u

    return [y1n, y2n]
end

Objective Function

The objective function consists of running cost and a fixed cost. The running cost function also has syntax similar to the dynamics function.

function L(x, u, t, p)
    return u[1]^2 + x[1]^2 + x[2]^2
end

The Final cost function involves the contribution of final state in the objective No final cost is involved

function phi(xf, uf, tf, p)
    return 0.0
end

Integral functions

Some of the problems can have integral constraints associated with them in each phase Since this problem does not have an integral constraint the integralfun will return nothing.

function integralfun(x, u, t, p)
    return nothing
end

Path functions

Some of the problems can have path constraints associated with them in each phase Since this problem does not have an path constraint the pathfun will return nothing.

function pathfun(x, u, t, p)
    y1 = x[1]; y2 = x[2]
    u = u[1]
    return -y2 + p.p1
end

function integralfun(x, u, t, p)
    return nothing
end

Now let us assign the various functions defined above to the phase ph that we have created

ph.L = L      # Adding the running cost
ph.phi = phi  # Adding the final time cost
ph.dyn = dyn  # Add the dynamics
ph.integralfun = integralfun # Add the integral constraint function
ph.pathfun = pathfun # Add the path constraint function
ph.n = n    # Number of points in initial mesh
ph.ns = ns  # State dimension
ph.nu = nu  # Input dimension
ph.nq = 0   # Dimension of the quadrature (integral) constraint
ph.np = 1   # Dimension of the path constraints
ph.p = p    # Auxillary parametrs named tuple

Let us select some of the options for the phase

Collocation method: Two options ["hermite-simpson", "trapezoidal"]. Default is "hermite-simpson"
Scale: Two options true or false

ph.collocation_method = "hermite-simpson"
ph.scale_flag = true

Now let us set the upper bounds and lower bounds on all the variables pathconstaints and integral constraints. The upper and lower bounds are defined in the limits field of the PH structure. The limits feild contains ll structure which corrsponds to lower limits. The ul corresponds to upper bounds. The ll and ul structures both contain the JuMP variables on which we wanrt to apply upper and lower bounds. These are: u : Input x : State xf : Final State xi: Initial State tf: Final State ti: Initial time If we want a particular variable to have a fixwd value set both the upper limits and the lower limits to the same value

ph.limits.ll.u = [-2.0]      # Lower bounds on input. Vector of dimension `nu`
ph.limits.ul.u = [2.0]  # Upper bounds on input. Vector of dimesion `nu`
ph.limits.ll.x = [-2.0, -2.0] # Lower bounds on state trajectory. Vector of dimension `ns`
ph.limits.ul.x = [2.0, 2.0] # Upper bounds on state trajectory. Vector of dimesion `ns`
ph.limits.ll.xf = [-2.0, -2.0]      # Lower bounds on final state. Vector of dimension `nu`
ph.limits.ul.xf = [2.0, 2.0]   # Upper bounds on input. Vector of dimesion `ns`
ph.limits.ll.xi = p.x0  # Lower bounds on initial state. Vector of dimension `ns`
ph.limits.ul.xi = p.x0  # Upper bounds on initial state. Vector of dimesion `ns`
ph.limits.ll.ti = 0.0   # Lower bounds on initial time. A Scalar.
ph.limits.ul.ti = 0.0   # Upper bounds on initial time. A Scalar
ph.limits.ll.tf = 5.0   # Lower bounds on input. Vector of dimension `nu`
ph.limits.ul.tf = 5.0   # Upper bounds on input. Vector of dimesion `ns`
ph.limits.ll.dt = 0.0     # Lower bounds on input. Vector of dimension `nu`
ph.limits.ul.dt = 5.0     # Upper bounds on input. Vector of dimesion `ns`
ph.limits.ll.k = [] # Lower bounds on input. Vector of dimension `nu`
ph.limits.ul.k = [] # Upper bounds on input. Vector of dimesion `ns`
ph.limits.ll.path = [-10.0]
ph.limits.ul.path = [0.0]

Set intial values

There are two options here "Auto" and "Manual". If "Auto" option is selected one need not specify tau and other init values

ph.set_initial_vals = "Auto"
ph.tau = range(start = 0, stop = 1, length = ph.n)
ph.xinit = ones(ph.ns, ph.n)
ph.uinit = ones(ph.nu, ph.n)
ph.tfinit = ph.limits.ll.tf
ph.tiinit = ph.limits.ll.ti
ph.kinit  = (ph.limits.ll.k + ph.limits.ul.k)/2

Specify global parameters

This problem only contains a single phase. However, in problems with multiple phases there are parameters which are global to all the phases. We specify it here.

Set limits on objective function

It has been observed that it is better to not set it as in keep the value false. However, an option is provided to set upper and lower values for the objective function

OC.set_obj_lim = false
OC.obj_llim = -2.0
OC.obj_ulim = 2.0

Setting up global parameters

OC.nkg: Number of global parameters
OC.kg_llim: Lower bound on global parameters
OC.kg_ulim: Upper bound on global parameters

Note that these have to be passed to the function through the auxillary parameters tupple p

OC.kg_llim = []
OC.kg_ulim = []

Setting up the event function

*OC.npsi: Number of constraints in event function *OC.psi_llim: Lower bound on constraint function *OC.psi_ulim: Upper bound on constraint function *OC.psi : Function which contains the event constraints

OC.npsi = 0
OC.psi_llim = []
OC.psi_ulim = []

Not the format of he event function psi. It takes the OCP object as input. The OCP object has all the phases stored in it in the field OC.ph. Each phase has state variables which can be accessed by ph.x and input variable which can be accessed by ph.u. In this problem we want that u must have a zero value at end of the phase.

function psi(ocp::DOC.OCP)
    (;ph) = ocp

    return nothing
end
OC.psi = psi

Call function to setup the JuMP model for solving optimal control problem

DOC.setup_mpocp(OC)

Solve for the control and state

DOC.solve_mpocp(OC)

DOC.solve(OC)

solution_summary(OC.model)

Compute co-states

lambda = DOC.differential_constraints_adjoints(ph, OC)
mu = DOC.path_constraint_adjoints(ph, OC)

Display results

println("Objective Value: ", objective_value(OC.model))

# using GLMakie
# f = Figure()
# ax1 = Axis(f[1,1])
# lines!(ax1, value.(ph.t), value.(ph.x[1,:]))
# ax2 = Axis(f[2,1])
# lines!(ax2, value.(ph.t), value.(ph.x[2,:]))
# # ax3 = Axis(f[1, 2])
# # lines!(ax3, value.(ph.t), value.(ph.x[3,:]))
# ax4 = Axis(f[2, 2])
# lines!(ax4,value.(ph.t), value.(ph.u[1,:]))
# display(f)

# fl = Figure()
# ax1 = Axis(fl[1,1])
# lines!(ax1, value.(ph.t[1:end-1]) ,lambda[1,:])
# ax2 = Axis(fl[2,1])
# lines!(ax2, value.(ph.t[1:end-1]) ,lambda[2,:])

# fmu = Figure()
# axmu = Axis(fmu[1,1])
# lines!(axmu, value.(ph.t), mu[1,:])

This page was generated using Literate.jl.