ControlSystem type is an interesting system because, unlike other system types, it cannot be numerically solved on its own. Instead, it must be transformed into another system before solving. Standard methods such as the "direct method", "multiple shooting", or "discretize-then-optimize" can all be phrased as symbolic transformations to a
ControlSystem: this is the strategy of this methodology.
Here we will start by defining a classic optimal control problem. Let:
where we want to optimize our controller
u(t) such that the following is minimized:
where $i$ is measured on (0,8) at 0.01 intervals. To do this, we rewrite the ODE in first order form:
is our loss function on the first order system.
Defining such a control system is similar to an
ODESystem, except we must also specify a control variable
u(t) and a loss function. Together, this problem looks as follows:
using ModelingToolkit @variables t x(t) v(t) u(t) @parameters p[1:2] D = Differential(t) loss = (4-x)^2 + 2v^2 + u^2 eqs = [ D(x) ~ v - p*x D(v) ~ p*u^3 + v ] sys = ControlSystem(loss,eqs,t,[x,v],[u],p)
One common way to solve nonlinear optimal control problems is by transforming them into an optimization problem by performing a Runge-Kutta discretization of the differential equation system and imposing equalities between variables in the same steps. This can be done via the
runge_kutta_discretize transformation on the
ControlSystem. While a tableau
tab can be specified, it defaults to a 5th order RadauIIA collocation, which is a common method in the field. To perform this discretization, we simply need to give a
dt and a timespan on which to discretize:
dt = 0.1 tspan = (0.0,1.0) sys = runge_kutta_discretize(sys,dt,tspan)
sys is an
OptimizationSystem which, when solved, gives the values of
u(t). Thus we solve the
OptimizationSystem using GalacticOptim.jl:
u0 = rand(length(states(sys))) # guess for the state values prob = OptimizationProblem(sys,u0,[0.1,0.1],grad=true) using GalacticOptim, Optim sol = solve(prob,BFGS())
And this is missing some nice interfaces and ignores the equality constraints right now so the tutorial is not complete.