Dynamic optimization

This web page is a duplicate of https://optimization.mccormick.northwestern.edu/index.php/Dynamic_optimization

Authors: Hanyu Shi (ChE 345 Spring 2014)

Steward: Dajun Yue, Fengqi You

Date Presented: Apr. 10, 2014

Authors: Issac Newton, Albert Einstein (ChE 345 Spring 2014)

Steward: Dajun Yue, Fengqi You

Date Presented: Apr. 10, 2014

Introduction

In this work, we will focus on the “at the same time” or direct transcription approach which allow a simultaneous method for the dynamic optimization problem. In particular, we formulate the dynamic optimization model with orthogonal collocation methods. These methods can also be regarded as a special class of implicit Runge–Kutta (IRK) methods. We apply the concepts and properties of IRK methods to the differential equations directly. With locating potential break points appropriately, this approach can model large-scale optimization formulations with the property of maintaining accurate state and control profiles. We mainly follows Biegler's work.

General Dynamic Optimization Problem

Differential algebraic equations in process engineering often have following characteristics: first,large-scale models – not easily scaled; second, sparse but no regular structure; third, direct linear solvers widely used; last, coarse-grained decomposition of linear algebra.

Figure 2. Dynamic optimization approach

There are several approaches can be applied to solve the dynamic optimization problems, which are shown in Figure 2.

Differential equations can usually be used to express conservation Laws, such as mass, energy, momentum. Algebraic equations can usually be used to express constitutive equations, equilibrium, such as physical properties, hydraulics, rate laws. Algebraic equations usually have semi-explicit form and assume to be index one i.e., algebraic variables can be solved uniquely by algebraic equations.

Dynamic Optimization Problem has the following general form:

${\begin{array}{l}\min \;\Phi \left({z\left(t\right),y\left(t\right),u\left(t\right),p,{t_{f}}}\right)\\s.t.\;\;{\frac {dz\left(t\right)}{dt}}=f\left({z\left(t\right),y\left(t\right),u\left(t\right),p}\right)\\g\left({z\left(t\right),y\left(t\right),u\left(t\right),p}\right)=0\\{z^{0}}=z\left(0\right)\\{z^{l}}\leq z\left(t\right)\leq {z^{u}}\\{y^{l}}\leq y\left(t\right)\leq {y^{u}}\\{u^{l}}\leq u\left(t\right)\leq {u^{u}}\\{p^{l}}\leq p\leq {p^{u}}\end{array}}$

$t$, time

$z$, differential variables y, algebraic variables

$t_{f}$ , final time

$u$, control variables

$p$, time independent parameters

(This follows Biegler's slides ）

Derivation of Collocation Methods

We first consider the differential algebraic system shown as follows:

${\begin{array}{l}{\frac {dz}{dt}}=f\left({z\left(t\right),y\left(t\right),u\left(t\right),p}\right),\;z\left(0\right)={z_{0}}\\g\left({z\left(t\right),y\left(t\right),u\left(t\right),p}\right)=0\end{array}}$ （1）

The simultaneous approach requires discretizing of the state variables $z\left(t\right)$, output variables $y\left(t\right)$ and manipulate variables $u\left(t\right)$. We require the following properties to yield an efficient NLP formulation:

1) The explicit ODE discretization holds little computational advantage because Since the nonlinear program requires an iterative solution of the KKT conditions.

2) A single step approach which is self-starting and does not rely on smooth profiles that extend over previous time steps is preferred, because the NLP formulation needs to deal with discontinuities in control profiles.

3) The high-order implicit discretization provides accurate profiles with relatively few finite elements. As a result, the number of finite elements need not be excessively large, particularly for problems with many states and controls.

Figure 1: Polynomial approximation for state profile across a finite element.

Polynomial Representation for ODE Solutions

We consider the following ODE:

${\frac {dz}{dt}}=f\left({z\left(t\right),t}\right),\;z\left(0\right)={z_{0}}$ (2)

to apply the collocation method, we need to solve the differential equation (2) at certain points. For the state variable, we consider a polynomial approximation of order $K+1$ (i.e., degree ≤ $K$ ) over a single finite element, as shown in the above figure. This polynomial, denoted by ${z^{K}}(t)$, can be represented in a number of equivalent ways, including the power series representation shown in equation (3), the Newton divided difference approximation, or B-splines.

${\frac {dz}{dt}}=f\left({z\left(t\right),t}\right),\;z\left(0\right)={z_{0}}$ (3)

We apply representations based on Lagrange interpolation polynomials to generate the NLP formulation, because the polynomial coefficients and the profiles have the same variable bounds. Here we select $K+1$ interpolation points in element i and represent the state in a given element $i$ as

${\begin{array}{l}\left.{\begin{array}{*{20}{c}}{t={t_{i-1}}+{h_{i}}\cdot \tau ,}\\{{z^{K}}\left(t\right)=\sum \limits _{j=0}^{K}{{l_{j}}\left(\tau \right)\cdot {z_{ij}},}}\end{array}}\right\}t\in \left[{{t_{i-1}},{t_{i}}}\right],\tau \in \left[{0,1}\right],\\where\;{l_{j}}\left(\tau \right)=\prod \limits _{k=0,\neq j}^{K}{\frac {\left({\tau -{\tau _{k}}}\right)}{\left({{\tau _{j}}-{\tau _{k}}}\right)}},\end{array}}$ （4）

${\tau _{0}}=0,\;{\tau _{i}}<{\tau _{i+1}},\;j=0,...,K-1$， and hi is the length of element $i$. This polynomial representation has the desirable property that ${z^{K}}({t_{ij}})={z_{ij}}$, where ${t_{ij}}={t_{i-1}}+{\tau _{j}}{h_{j}}$.

We use a Lagrange polynomial with K interpolation points to represent the time derivative of the state. This leads to the Runge–Kutta basis representation for the differential state:

${z^{K}}\left(t\right)={z_{i-1}}+{h_{i}}\cdot \sum \limits _{j=1}^{K}{{\Omega _{j}}\left(\tau \right)}\cdot {{\dot {z}}_{ij}}$ (5)

where ${z_{i-1}}$ is a coefficient that represents the differential state at the beginning of element $i$, ${{\dot {z}}_{ij}}$represents the time derivative ${\frac {d{z^{K}}({t_{ij}})}{d\tau }}$, and ${\Omega _{j}}(\tau )$ is a polynomial of order K satisfying

${\Omega _{j}}\left(\tau \right)=\int _{0}^{\tau }{{l_{j}}(\tau ')}d\tau ',t\in \left[{{t_{i-1}},{t_{i}}}\right],\tau \in \left[{0,1}\right]$ (6)

We substitute the polynomial into equation (1) to calculate the polynomial coefficients, which is an approximation of the DAE. This results in the following collocation equations.

${z^{K}}\left(t\right)={z_{i-1}}+{h_{i}}\cdot \sum \limits _{j=1}^{K}{{\Omega _{j}}\left(\tau \right)}\cdot {{\dot {z}}_{ij}}$ (7)

with ${z^{k}}({t_{i}}-1)$ calculated separately. For the polynomial representations (4) and (5), we normalize time over the element, write the state profile as a function of τ , and apply ${\frac {d{z^{K}}}{d\tau }}={h_{i}}{\frac {d{z^{K}}}{dt}}$ easily. For the Lagrange polynomial (4), the collocation equations become

$\sum \limits _{j=0}^{K}{{z_{ij}}\cdot {\frac {d{l_{j}}\left({\tau _{k}}\right)}{d\tau }}}={h_{i}}\cdot f\left({{z_{ik}},{t_{ik}}}\right),\;k=1,...,K$ (8)

while the collocation equations for the Runge–Kutta basis are given by

${{\dot {z}}_{ik}}=f\left({{z_{ik}},{t_{ik}}}\right)$ (9)

${z_{ik}}={z_{i-1}}+{h_{i}}\cdot \sum \limits _{j=1}^{K}{{\Omega _{j}}\left(\tau \right)}\cdot {{\dot {z}}_{ij}},\;k=1,...,K$ (10)

with ${z_{i}}-1$determined from the previous element $i-1$ or from the initial condition on the ODE.

Example

An example is given here to demonstrate the application of the collocation method.

A differential equation is given as follows:

${\frac {dz}{dt}}={z^{2}}-2\cdot z+1,\;z\left(0\right)=-3$ (11)

With t \in \left[ {0,1} \right], The analytic solution of this differential equation is $z\left(t\right)={\frac {4\cdot t-3}{4\cdot t+1}}$.

Lagrange interpolation and collocation method is applied to this differential equation respectively. And the number of collocation points in each finite element is 3. The number of finite elements is N , and the length of the finite element is 1/N . The following equations is given then:

$\sum \limits _{j=0}^{3}{{z_{ij}}{\frac {d{l_{j}}\left({\tau _{k}}\right)}{d\tau }}}=h\left({z_{ik}^{2}-2\cdot {z_{ik}}+1}\right),\;k=1,...,3,\;i=1,...,N$ (12)

${z_{i+1,0}}=\sum \limits _{j=0}^{0}{{l_{j}}\left(1\right)}\cdot {z_{ij}},\;i=1,...,N-1$ (13)

${z_{f}}=\sum \limits _{j=0}^{K}{{l_{j}}\left(1\right)}\cdot {z_{Nj}},\;{z_{1,0}}=-3$ (14)

With Radau collocation method, ${\tau _{0}}=0$, ${\tau _{1}}={\rm {0.155051}}$, ${\tau _{2}}={\rm {0.644949}}$ and ${\tau _{3}}=1$ can be obtained. The collocation equations are given as follows:

$\sum \limits _{j=0}^{3}{{z_{j}}{\frac {d{l_{j}}\left({\tau _{k}}\right)}{d\tau }}}=\left({z_{k}^{2}-2\cdot {z_{k}}+1}\right),\;k=1,...,3$ (14)

which can be formulated as:

${\begin{array}{l}{z_{0}}\cdot \left({-30\cdot \tau _{k}^{2}+36\cdot {\tau _{k}}-9}\right)+{z_{1}}\cdot \left({{\rm {46.7423}}\cdot \tau _{k}^{2}-{\rm {51.2592}}\cdot {\tau _{k}}+{\rm {10.0488}}}\right)\\+{z_{3}}\cdot \left({-{\rm {26.7423}}\cdot \tau _{k}^{2}+{\rm {20.5925}}\cdot {\tau _{k}}-{\rm {1.38214}}}\right)+{z_{3}}\cdot \left({10\cdot \tau _{k}^{2}-{\frac {16}{3}}\cdot {\tau _{k}}+{\frac {1}{3}}}\right)\\=\left({z_{k}^{2}-2\cdot {z_{k}}+1}\right),\;k=1,...,3\end{array}}$ (15)

Figure 3. Comparison of Radau collocation solution with exact solution

The results are given as following by solving the above equations:

$\left\{{\begin{array}{*{20}{c}}{{z_{0}}=-3}\\{{z_{1}}=-{\rm {1.65701}}}\\{{z_{2}}={\rm {0.032053}}}\\{{z_{3}}={\rm {0.207272}}}\end{array}}\right.$ (16)

As shown in Figure 3.2 the error $\left\|{z\left(1\right)-{z^{K}}\left(1\right)}\right\|$， is less than $10^{-}6$ for $N=5$ and converges with $O(h^{5})$, which is consistent with the expected order $2K-1$.

(This example follows the work of Biegler and can be found in P293 of “ Nonlinear Programmng”.)

Conclusion

In this work, we mainly discussed simultaneous collocation approach for dynamic optimization problems, which formulated the differential equations to a set of algebraic equations. These direct transcription formulations depended on fully discretizing of the differential algebraic equations (DAE), which enabled us solve the simultaneous optimization problem without relying on embedded DAE solvers. Because of this simultaneous formulation, we got the exact first and second order derivatives through the optimization modeling system, and both structure and sparsity can be exploited.

References

1. Biegler, Lorenz T. Nonlinear programming: concepts, algorithms, and applications to chemical processes. Vol. 10. SIAM, 2010.

2. Chu, Yunfei, and Fengqi You. "Integration of scheduling and control with online closed-loop implementation: Fast computational strategy and large-scale global optimization algorithm." Computers & Chemical Engineering 47 (2012): 248-268.

3. http://en.wikipedia.org/wiki/Dynamic_programming

4. http://en.wikipedia.org/wiki/Differential_algebraic_equation

5. http://numero.cheme.cmu.edu/uploads/dynopt.pdf