Optimization¶
We now describe the interface for defining and solving switching time optimization problems.
Problem Definition¶
This package allows us to define and solve problems in the form
where the decision variable is the vector of \(N+1\) intervals \(\delta = \begin{bmatrix}\delta_0 & \dots & \delta_{N}\end{bmatrix}^\top\in \mathbb{R}^{N+1}\) such that \(\delta_i = \tau_{i+1} - \tau_i\). Each interval \(\delta_i\) defines how long the \(i\)-th dynamics are active. The state trajectory is \(x(t) \in \mathbb{R}^{n}\). The value \(T_\delta\) is defined as the final time when intervals \(\delta\) are applied, i.e.
The set \(\Delta\) defines the set of feasible intervals
The variable \(T\) defines the desired final time of the interval. The scalars \(lb_i\) and \(ub_i\) define additional constraints on the interval in case we would like to have a minimum or a maximum time in which the \(i\)-th dynamics are active.
Linear Dynamics¶
In case when the dynamics are linear of the form
we can define a 3-dimensional matrix A
whose slices A[:,:,i]
represent dynamics \(A_{i-1}\):
A = Array(Float64, n, n, N+1)
A[:, :, i] = ... # Dynamics A_{i-1}
...
The switching time optimization problem can be quickle defined as
p = stoproblem(x0, A)
Where x0
is the initial state vector \(x_0\) and A
is the 3-dimensional matrix defining the dynamics.
Noninear Dynamics¶
Given a nonlinear system defined by dynamics
where \(u(t)\) the input vector assuming integer values \(u_i\) between switching instants
we can define our switched nonlinear system as
To create the optimization problem we need to define the nonlinear dynamics by means of an additional function
function nldyn(x, ui)
...
end
returning the vector of states derivatives. The variable x
is the state \(x(t)\) and ui
is the input vector \(u_i\). Moreover, we need to define the jacobian of the switched dynamics with respect to the system states
by means of an additional function
function jac_nldyn(x, ui)
...
end
Note that the function jac_nldyn
returns a matrix having in each row the gradient of every component of the function \(f_i(x(t))\) with respect to each state component. Last element necessary to construct the matrix U
having a column each integer input vector ui
. Then, we can define the switching time optimization problem as:
p = stoproblem(x0, nldyn, jac_nldyn, U)
Note
The nonliner switched system optimization operates by introducing additional linearization points at an equally spaced linearization grid. To set the number of linearization points to \(100\) for example, it is just necessary to add an extra argument to the previous function call as follows:
p = stoproblem(x0, nldyn, jac_nldyn, U, ngrid = 100)
where ngrid
defines the number of linearization points.
Optional Arguments¶
There are many additional keyword arguments that can be be passed to the stoproblem(...)
function to customize the optimization problem.
Parameter | Description | Default value |
---|---|---|
t0 |
Initial Time \(t_0\) | 0.0 |
tf |
Final Time \(t_0\) | 1.0 |
Q |
Cost matrix \(Q\) | eye(n) |
lb |
Vector of lower bounds \(lb_i\) | zeros(N+1) |
ub |
Vector of lower bounds \(ub_i\) | Inf*ones(N+1) |
tau0ws |
Warm starting initial switching times | Equally spaced between t0 and tf |
solver |
MathProgbase.jl solver | IpoptSolver() |
Problem Solution¶
Once the problem is defined, it can be solved by simply running
solve!(p)
Choosing Solver¶
Any NLP solver supported by JuliaOpt may be used through MathProgBase.jl interface. The default solver is Ipopt. To use KNITRO solver with the linear example, it is just necessary to specify an AbstractMathProgSolver
object (see here for more details) when the problem is created
using KNITRO
p = stoproblem(x0, A, solver = KnitroSolver())
All the solver-specific options can be passed when creating the AbstractMathProgSolver
object: algorithm types (first/second order methods), tolerances, verbosity and so on.
Obtaining Results¶
The optimal cost function and the optimal switching times and intervals can be obtained as follows:
objval = getobjval(p)
tauopt = gettau(p)
deltaopt = getdelta(p)
We can get the execution time (including the time for the function calls) and the status of the solver by executing:
stat = getstat(p)
soltime = getsoltime(p)
Optimizing in a Loop¶
The toolbox is suited for receeding horizon implementations. To run the optimization in a loop it is just necessary to update the value of the current state x0
and to update the warm starting point tau0ws
which is usually chosen as the optimal solution at the previous optimizaton.
To set the initial state at x0
it is just necessary to return
setx0!(m, x0)
We can set the warm starting point at tau0ws
with
setwarmstart!(m, tau0ws)