NUMERICAL OPTIMAL CONTROL OF A COUPLED ODE-PDE MODEL OF A TRUCK WITH A FLUID BASIN

. We consider a numerical study of an optimal control problem for a truck with a ﬂuid basin, which leads to an optimal control problem with a coupled system of partial diﬀerential equations (PDEs) and ordinary diﬀerential equations (ODEs). The motion of the ﬂuid in the basin is modeled by the nonlinear hyperbolic Saint-Venant (shallow water) equations while the vehicle dynamics are described by the equations of motion of a mechanical multi-body system. These equations are fully coupled through boundary conditions and force terms. We pursue a ﬁrst-discretize-then-optimize approach using a Lax-Friedrich scheme. To this end a reduced optimization problem is obtained by a direct shooting approach and solved by a sequential quadratic programming method. For the computation of gradients we employ an eﬃcient adjoint scheme. Numerical case studies for optimal braking maneuvers of the truck and the basin ﬁlled with a ﬂuid are presented.

Abstract. We consider a numerical study of an optimal control problem for a truck with a fluid basin, which leads to an optimal control problem with a coupled system of partial differential equations (PDEs) and ordinary differential equations (ODEs). The motion of the fluid in the basin is modeled by the nonlinear hyperbolic Saint-Venant (shallow water) equations while the vehicle dynamics are described by the equations of motion of a mechanical multi-body system. These equations are fully coupled through boundary conditions and force terms. We pursue a first-discretize-then-optimize approach using a Lax-Friedrich scheme. To this end a reduced optimization problem is obtained by a direct shooting approach and solved by a sequential quadratic programming method. For the computation of gradients we employ an efficient adjoint scheme. Numerical case studies for optimal braking maneuvers of the truck and the basin filled with a fluid are presented. reference system, d W the distance of the basin relative to the (d, z) reference system, andḋ T andḋ W its velocities. A vertical motion of the truck does not apply in this model. Assuming linear spring and damper force laws, the spring-damper force F , which acts between the truck and the basin, is given by where c is the spring force constant, k is the damper force constant, andd is an offset used in the spring force. Ifd = d W (0) − d T (0), then the basin is at rest at time t = 0 relative to the truck. The dynamics of the system are given by the following coupled system of ordinary differential equations (ODEs) and partial differential equations (PDEs), wherein the ODE is given by Herein, friction terms are not considered for simplicity. The force term m W a W is due to the average acceleration of the basin owing to fluid motion. It is defined by where v t (t, x) denotes the acceleration of the fluid at time t at position x. The motion of the fluid in the basin is given by the 1d Saint-Venant equations with Ω := (0, T ) × (0, L) and initial and boundary conditions The equations (4) and (5) prescribe conservation of mass and conservation of (linear) momentum. Under the assumption that typical vertical scales are much smaller than horizontal scales, allowing to average over the depth of the fluid, the Saint-Venant equations may be derived [11,Ch. 2] from the Navier-Stokes equations with zero normal flow and shear stresses at boundaries. Contrary to the case of the Boussinesq equations, no frequency dispersion of the waves is modeled here. Note that the fluid level at (t, x) is given by B(x) + h(t, x). For simplicity, we assume that 0 < h ≤ h(t, x) ≤ h < H(x) holds for all (t, x) ∈Ω, H(x) describing the given wall height (or the cap) of the basin.
Subject to the assumption h = 0, equation (5) can be simplified to the so-called nonconservative form We remark that a velocity is not subject to a fundamental conservation principle. Therefore (6) does not hold across jumps in general. Boundary conditions for h at x = 0 and x = L can be obtained from (6), if it is considered at The average acceleration (3) of the basin computes to where the boundary conditions v(t, 0) = v(t, L) = 0 have been exploited. With this, the spring-damper force in Eq.
(2) cancels out and it simplifies to 2. An optimal control problem for the coupled ODE-PDE system. We consider optimal control problems with free final time T > 0 for the motion of the truck from a given initial state to a given terminal state defined by the boundary conditions d W are given numbers. A collection of potential objective functions J i , i = 0, 1, 2, 3, is given as follows: • final time J 0 := T to be minimized; • minimization of the deviation of a given fluid level h d (typically constant): • minimization of the control effort: • minimization of the distance of terminal velocities to given terminal velocities of the truck and the basin: Combining these objectives with suitable non-negative constants α i ≥ 0, i = 0, . . . , 3, yields the objective function J := α 0 J 0 + α 1 J 1 + α 2 J 2 + α 3 J 3 to be minimized subject to the dynamics (1), (2), (4), (6) and the initial and boundary conditions. In addition, bounds on the control u and the state h can be added, i.e.
For its numerical solution, it is convenient to transform the optimal control problem with free final time T > 0 to an equivalent optimal control problem on a fixed time horizon, e.g. on [0, 1]. This is achieved by the linear time transformation t(τ ) := τ T with τ ∈ [0, 1] and T > 0 being an additional optimization variable. Moreover, we write the ODE as a first order system by introducing the velocities η T :=ḋ T and η W :=ḋ W . Application of these transformations yields the following optimal control problem (OCP) on the fixed time interval [0, 1]. For notational convenience the time is denoted by t again and not by τ . Likewise, the transformed states are denoted by

OCP: Minimize
withΩ := (0, 1) × (0, L) and the initial and boundary conditions A related control problem (without a truck) with the aim to stabilize a given water level profile h d was considered in [1, Sect. 6.3].
3. Discretization of the Saint-Venant equations and the ODEs. We apply a direct discretization technique to the optimal control problem OCP in Section 2. Herein, the Saint-Venant equations are being discretized by the Lax-Friedrich (LF) scheme, which is an explicit method of first order in time and of second order in space. As a consequence of the explicit scheme, the Courant-Friedrichs-Levy (CFL) condition has to be obeyed in order to guarantee stability of the discretized solution. Moreover, the CFL condition is necessary and sufficient for the convergence of the LF scheme [10,Sect. 8.3]. The LF scheme is appealing owing to its simplicity (in contrast to implicit schemes), the availability of convergence results (in contrast to, e.g., ENO schemes), and because it introduces artificial diffusion to the problem, which matches well with the discussion of viscosity solutions.
We discretize the closure of the domainΩ by an equidistant grid in time and space with grid points assuming that the CFL condition is satisfied. It turns out numerically that seems to be an appropriate choice, see Section 3.3 below. Define the approximations for the function q := 1 2 v 2 + gh. Application to the Saint-Venant equations yields with which together with the boundary conditions (7) and (8) yields Discretization of the ODE states d T , d W , η T , and η W by the explicit Euler method on the time grid yields the difference equations It is convenient to write these discretization schemes in the compact form z n+1 = z n + ∆t Φ(t n , z n , u n , T ), n = 0, 1, . . . , N − 1, where z n := (d n T , d n W , η n T , η n W , h n 1 , . . . , h n M −1 , v n 1 , . . . , v n M −1 ) for n = 0, 1, . . . , N denotes the discrete state at time point t n and Φ is the increment function, which is defined by the LF scheme for the PDE part and the explicit Euler scheme for the ODE part. Solving the recursion (12) step-by-step for z n yields the relation z n = z n (p) := (d n T (p), d n W (p), η n T (p), η n W (p), h n 1 (p), . . . , h n M −1 (p), v n 1 (p), . . . , v n M −1 (p)) with p := (u 0 , u 1 , . . . , u N −1 , T ) , i.e. the discrete state trajectory depends only on the discrete control input u 0 , . . . , u N −1 and the final time T (single shooting idea). Notice, that the initial state is fixed. Approximation of the integrals in the objective function J of OCP by the second order trapezoidal rule in the space variable x and by first order Riemann sums in the time variable t yields the following reduced discretized optimal control problem (RDOCP): RDOCP: Minimizẽ The finite dimensional optimization problem RDOCP can be solved numerically by, e.g., the software package OCPID-DAE1, see [5]. It uses a gradient based sequential quadratic programming method with BFGS Hessian approximations and thus, gradients of the objective function and the nonlinear constraints have to be provided in an efficient way. To this end we employ an adjoint method. This method is particularly efficient if only few nonlinear constraints are present, for instance if the state constraints h ≤ h n i (p) ≤ h could be neglected. We explain merely the basic principle of the adjoint method and leave explicit derivations to the reader, see [4, Section 5.3] for details. We intend to compute the gradient of some function of type Γ(p) := γ(z n (p), p), where p = (u 0 , . . . , u N −1 , T ) and n ≤ N . Notice that the functions d N T (p), d N W (p) in the terminal constraints, and the functions h n i (p) in the state constraints fit into this form, since those are components of z N and z n , respectively. The objective function does not yet fit into this form, but it can be transformed into a Mayer-type objective function by adding an artificial state and a corresponding discrete state equation, which fits into the form (12). This transformation will allow to apply the following adjoint approach to the modified objective function and to compute its gradient. Using a formal Lagrange technique one can derive the expression with the increment function Φ from (12), compare [4, Section 5.3]. Herein, λ solves the discrete adjoint equation Notice that the adjoint equation is solved backwards in time. Moreover, the adjoint equation has to be re-solved for every function for which a gradient has to be computed. If no state constraints are present in RDOCP, then the gradients of the objective function and the two boundary conditions need to be computed, i.e. only three adjoint equations need to be solved to obtain the desired gradients independently on the dimension of p. This independence of the dimension of p is particularly nice in our case since the CFL condition (10) usually requires to use many time grid points, which leads to a high number of controls and thus a high dimension of p. If the state constraints are present in RDOCP, then the adjoint approach is not very efficient anymore, since an adjoint equation needs to be solved for each of the N(M+1) state constraints (not all on the full time horizon, though). In this case a standard sensitivity analysis approach for the discretized dynamics or even finite differences are more efficient. with an artificial viscosity µ = ∆x 2 2T ∆t , which is introduced by the Lax-Friedrich scheme (analogously for the v-equation). The CFL condition (10) requires ∆t and ∆x to be chosen such that holds [10, p. 192 (10) has turned out to be appropriate and yields µ = 30 L 2T ∆x. Hence, for small ∆x, the modification of the original equation is small. The Laplacian has a regularization effect and avoids numerical difficulties, if the initial data h 0 is sufficiently smooth [8,Example 2.2.6]. In general, hyperbolic conservation laws exhibit shocks or rarefaction waves even for smooth initial data [3,Sect. 11.3.2]. Finite differences are a special case of finite volumes and preserve conservation laws in the discretized version. The LF scheme exhibits consistence of the numerical flux with the conservation law and it is consistent with a so-called entropy condition. Thus, according to the Lax-Wendroff theorem [8,Th. 4.4.1] applied to a hyperbolic system 1d in space, a bounded numerical solution of the LF scheme converges globally to a weak solution as ∆x → 0. Provided (h 0 , v 0 ) ∈ L ∞ ∩ L 1 (R), the LF scheme approximates (w.r.t. to weak- * convergence) the entropy solution (h, v) ∈ L ∞ (R) of a pure initial value problem [8, Example 2.2.6/Th. 3.3.28 ]. The convergence to an entropy solution is due to the numerical damping term. Other explicit schemes where the Lax-Wendroff theorem can be applied are e.g. the Glimm scheme or the Godunov scheme [8,Lemma 4.4.5]. There are many other possibilities to discretize the differential equations. We could use more refined semi-explicit schemes, e.g. MUSCL/Kurganov-Tadmor or Kurganov-Levy [9], or implicit schemes, e.g. ENO-II. Implicit schemes might allow to avoid the restrictive CFL condition that causes comparatively large computing times for the optimal control. For a comparison of explicit and implicit schemes, claiming that explicit schemes are competitive, see [8,Remark 2.3.34]. For the Kurganov-Levy scheme, it can be proved that h remains non-negative in the discretization.

3.3.
Courant-Friedrichs-Levy (CFL) condition. We justify theoretically our numerical estimate for the CFL constant (10). As in [2, Sect. 2.1] we rewrite the system (4) and (6) by diagonalization of the Jacobian J F for F = hv, 1 2 v 2 + gh as where R 1/2 := v ± 2 √ gh are the so-called Riemann invariants and σ 1/2 = v ∓ √ gh are the eigenvalues of the Jacobian J F . The eigenvalues determine the maximal group velocity V := max x,t max i=1,2 |σ i (x, t)|, i.e. the maximal speed of propagation. Furthermore, since the two eigenvalues are distinct, we see that our system (4) and (6)  The control bounds are u min = −20000 and u max = 2000, reflecting that the truck can brake faster than accelerate. We choose N = 1500, M = 50, the optimality tolerance 5 · 10 −7 , and the feasibility tolerance 10 −8 within the SQP method, as described in [5], for the first two test examples.
For an optimal emergency braking maneuver, we consider the minimum time problem with the weights α 0 = 0.1, α 1 = α 2 = 0 and α 3 = 200 and without state constraints for h. The numerical solution is depicted in Figure 2. The optimal objective function value

5.
Extensions, related work, and future research. Controllability of a similar system, but without the boundary conditions (7) and (8) and zero source terms B and F , has been investigated in [2,Sect. 3]). It turns out that the linearized system is not controllable in general, but steady-state controllable. The nonlinear system is locally controllable around the equilibrium point [1,Sect. 6.3]. Global boundary controllability of the Saint Venant equations between steady states is established in [6] and more generally for sloped canals with friction in [7].
A theoretical investigation of the optimal control problem OCP with regard to existence of optimal solutions and necessary conditions of optimality will be the focus of an upcoming study in order to underpin the numerical results of this paper. Likewise a parametric     sensitivity analysis based on [4, Section 6.1] allows to study the dependence of an optimal solution on model parameters (e.g. masses m T , m W , spring-damper coefficients c, k) or disturbances.
Moreover, the problem can be extended to a 3d problem with 2d Saint-Venant equations and a full 3d model of a truck. Instead of using a direct shooting approach it is also possible to use direct collocation and to exploit the sparse structure using exact derivatives. Preliminary numerical results are promising and even show a local quadratic convergence.