**4. Numerical solution of the optimal control problem**

The optimal control problem is numerically solved as large scale structured non-linear pro‐ gramming problem (NLP):

$$\min\_{\mathbf{y}} \quad \left\{ J(\mathbf{y}) \mid \mathbf{h}(\mathbf{y}) = \mathbf{0} ; \mathbf{g}(\mathbf{y}) \le \mathbf{0} \right\} \tag{15}$$

The optimization variables are the state and control variables of the several stages in time:

$$\mathbf{y} = \left[ \left( \mathbf{x}^{0} \right)^{\mathrm{T}} \left( \mathbf{u}^{0} \right)^{\mathrm{T}} \dots \left( \mathbf{x}^{K-1} \right)^{\mathrm{T}} \left( \mathbf{u}^{K-1} \right)^{\mathrm{T}} \left( \mathbf{x}^{K} \right)^{\mathrm{T}} \right]^{\mathrm{T}} \tag{16}$$

The state equations of the process model are incorporated as equality constraints:

$$\mathbf{h}(\mathbf{y}) = \begin{bmatrix} \mathbf{x}^0 - \mathbf{x}(t\_0) \\ \mathbf{f}^0(\mathbf{x}^0, \mathbf{u}^0, \mathbf{z}^0) - \mathbf{x}^1 \\ \vdots \\ \mathbf{f}^{K-1}(\mathbf{x}^{K-1}, \mathbf{u}^{K-1}, \mathbf{z}^{K-1}) - \mathbf{x}^K \end{bmatrix} \tag{17}$$

The advantage of this problem formulation with *(K(n+m)+n)* optimization variables is the special sparsity structure with a block-diagonal Hessian-matrix and block-banded Jacobian matrices (*k*: number of time steps, *n*: number of state variables, *m*: number of control varia‐ bles), which follows from the fact, that the process equations as coupling element between adjusting stages are linear in the state variables *xk+1*. The numerical solution of this problem requires about (*K(n+m)3* ) basic arithmetic operations.

One alternative approach consists of eliminating the state variables from the optimization problem. The reduced dimension of the according non-linear programming problem with (*K m*) optimization variables comes along with a loss of structure. The Hessian and Jacobian matrices are full. Because of the solution effort of order (*K m*)*<sup>3</sup>* and the large amount of con‐ trol variables (number of edges in the network description) this approach is not promising for this special field of application.

,min ,max *k kk*

,min ( ) ,max *k kk k*

With respect to the practical applicability selected parts of the inequality constraints (10) can be relaxed in order to avoid infeasible optimization problem with respect to unrealistic man‐

The optimal control problem is numerically solved as large scale structured non-linear pro‐

min ; {*J*( ) ( ) = £ ( ) } **<sup>y</sup>**

The optimization variables are the state and control variables of the several stages in time:

( ) ( ) ( ) ( ) ( ) TT T T 00 1 1 *<sup>T</sup> <sup>T</sup>* é ù *K KK* - - <sup>=</sup> ê ú

( )

*t*

0 0 0 00 1

( )

, ,

é ù - ê ú

0

**x x fxuz x**


) basic arithmetic operations.

( )

**fxuz x**


The advantage of this problem formulation with *(K(n+m)+n)* optimization variables is the special sparsity structure with a block-diagonal Hessian-matrix and block-banded Jacobian matrices (*k*: number of time steps, *n*: number of state variables, *m*: number of control varia‐ bles), which follows from the fact, that the process equations as coupling element between adjusting stages are linear in the state variables *xk+1*. The numerical solution of this problem

, , *K K KK K*

1 1 11


The state equations of the process model are incorporated as equality constraints:

( )

**h y**

and the hydraulic head *hhydr* of the observation wells:

40 Water Supply System Analysis - Selected Topics

**4. Numerical solution of the optimal control problem**

agement demands.

gramming problem (NLP):

requires about (*K(n+m)3*

*v vv* **x xx** £ £ (13)

*hydr gw hydr* **h xh** £ £ *g* (14)

**y hy 0gy 0** (15)

<sup>M</sup> (17)

ë û **yx u x u x** <sup>K</sup> (16)

For the numerical solution of large scale non-linear programming problems interior point (IP) solver have become popular during the last years because of their superior behavior for NLPs with many inequality constraints. In this approach the objective function is expanded by adding barrier terms for the inequality constraints:

$$\min\_{\mathbf{y}} \quad \left\{ J(\mathbf{y}) + \mu \sum\_{j=1}^{n\_{\mathbf{g}}} \ln \left( -\mathbf{g}\_{j}(\mathbf{y}) \right) \, \Big|\, \mathbf{h}(\mathbf{y}) = \mathbf{0} \right\} \tag{18}$$

The solution of the original NLP (15) results from the subsequent solution of (17) with a de‐ caying sequence of *µ → 0*. The identification of right active set with its combinatorial com‐ plexity is avoided. The state of the art non-linear interior point solver IPOPT is used for the application at hand [15]. The interface for multistage optimal control problems of the opti‐ mization solver HQP [16], which provides an efficient way for problem formulation along with routines for the derivation of ∇ *J*, ∇*h* , ∇*g* by means of automatic differentiation [17], is used and coupled to IPOPT.

A typical water management problem (horizon of five years, discretization of one month) has about 8000 optimization variables, 5500 equality constraints and 7200 inequality con‐ straints. The numerical solution takes approximately 60 iterations and a calculation time of 10 seconds on an Intel Core 2 Duo CPU (2.5 GHz). Table 1 shows the linear depend‐ ency of the numerical solution effort from the number of time steps within the optimiza‐ tion horizon.


**Table 1.** Numerical solution effort in dependence on the optimization horizon
