Optimization techniques for tree-structured nonlinear problems

Robust model predictive control approaches and other applications lead to nonlinear optimization problems defined on (scenario) trees. We present structure-preserving Quasi-Newton update formulas as well as structured inertia correction techniques that allow to solve these problems by interior-point methods with specialized KKT solvers for tree-structured optimization problems. The same type of KKT solvers could be used in active-set based SQP methods. The viability of our approach is demonstrated by two robust control problems.


Introduction
This paper addresses nonlinear optimization problems (NLPs) with an underlying tree topology.The prototypical example are multistage stochastic optimization problems.These are computationally expensive because they involve some random process whose discretization yields a scenario tree that grows exponentially with the length of the planning horizon.Primal and dual decomposition methods are specially tailored to linear and mildly nonlinear convex multistage stochastic optimization problems.Moreover, they extend well to the mixed-integer case.See Birge and Louveaux (2011), Kall and Wallace (1994) for an overview of algorithmic approaches and of stochastic optimization in general.To tackle highly nonlinear convex and nonconvex NLPs, we consider interior-point and SQP methods.Interior-point methods (IPMs) are also efficient on linear or quadratic problems, while SQP methods extend well to mixedinteger problems.In both cases the key to efficiency is the algebraic structure of the arising KKT systems rather than the stochastic structure.Therefore, we allow arbitrary trees instead of only scenario trees.All methods mentioned above are amenable to parallelization due to the rich structure of the underlying tree.
In Steinbach (2002) and the references therein, the third author has developed suitable formulations and a sequential interior-point approach for convex tree-structured NLPs with polyhedral constraints, although with parallelization and with the fully nonlinear and nonconvex case in mind.Other early interior-point approaches for stochastic optimization include Berger et al. (1995), Carpenter et al. (1993), Czyzyk et al. (1995), Jessup et al. (1994), Schweitzer (1998).For more recent structure-exploiting parallel algorithms see, e.g., Blomvall (2003), Blomvall and Lindberg (2002), Gondzio and Grothey (2006, 2007, 2009), Lubin et al. (2012).Based on the dissertation (Hübner 2016), a massively distributed implementation of our algorithm has been presented in Hübner et al. (2017).Here we address the extension to fully nonlinear and nonconvex problems.In particular, we present proper NLP formulations, structure-preserving Quasi-Newton updates, and tailored inertia correction techniques.Although all presented techniques possess a straightforward distributed implementation, we will not address this in detail.
The paper is structured as follows.In Sect. 2 we briefly introduce our approach and extend the problem formulations and KKT systems of Steinbach (2002) to the general NLP case.Section 3 presents (partially) separable Quasi-Newton update formulas that preserve the problem-specific block structure.The proper specialization of standard inertia correction techniques to this structure is discussed in Sect. 4. In Sect.5, finally, two case studies from robust process control that lead to multistage stochastic NLPs with ODE dynamics serve as proof of concept for our algorithmic techniques.

Tree-sparse optimization
In this paper, we consider the integrated modeling and solution framework for treestructured convex NLPs from Steinbach (2002) that we refer to as "tree-sparse".It consists of natural model formulations that have favorable regularity properties and block-level sparsity admitting O(|V |) KKT solution algorithms, where |V | is the number of tree nodes.There are three standard forms that cover, respectively, the general case with no further structure and two more specialized control formulations with different stochastic interpretations.The NLPs that we study below have linearizations whose polyhedral constraints match precisely the structure considered in Steinbach (2002).
For the following, let V = {0, . . ., N } denote the node set of a tree rooted in 0. Given a node j ∈ V , we denote the set of successors by S( j), the predecessor by i = π( j) (if j = 0), and the path to the root by Π( j) = ( j, π( j), . . ., 0).The level of j is the path length, i.e., t( j) = |Π( j)| − 1.Finally, L t stands for the set of level t nodes and L for the set of leaves.In a scenario tree we have L = L T where T is the tree's depth, and every node j has a probability p j such that j∈L t p j = 1 for every t.

Tree-sparse NLPs
Consider a general NLP with equality constraints, range inequality constraints, and simple bounds in the form (1) We call the NLP tree-sparse if it satisfies certain separability properties.Given a tree with vertex set V and variable vector y = (y j ) j∈V , the objective φ and range constraints c R have to be separable with respect to the node variables y j , and the equality constraints c E split into dynamic constraints with a Markov structure (y j depends only on its predecessor y i ) and separable global constraints.We distinguish three variants of tree-sparse NLPs depending on the type of dynamic constraints: implicit tree-sparse NLPs have implicit dynamics while outgoing and incoming tree-sparse NLPs with y j = (x j , u j ) and y j = (u j , x j ), respectively, have explicit dynamics in control form where the current state x j depends on (x i , u i ) or (x i , u j ).In the following, we state these three NLP variants without further explanation.All details can be found in Steinbach (2002).
The implicit tree-sparse NLP reads min the outgoing tree-sparse NLP is given by min x,u j∈V and the incoming tree-sparse NLP reads min u,x j∈V At the root, j = 0, the preceding variable y i is empty and the dynamic constraints reduce to respective initial conditions g 0 (y 0 ) = h 0 , x 0 = h 0 , and x 0 = h 0 (u 0 ).
In addition to global equality constraints one might also consider local equality constraints as in Steinbach (2002).In the above NLPs these would take the respective forms To avoid unnecessary technical complexity we assume that local constraints are modeled as global constraints.Specific discussions will be added where the distinction is relevant.Again for simplicity we consider range constraints in the form c R (y) ≥ 0 rather than lower and upper constraints as in Steinbach (2002), c R (y) ∈ [r − , r + ].

Tree-sparse KKT systems
We solve the NLP (1) by an interior-point method, handle bounds directly, and convert range constraints to c R (y) − s = 0 using slacks s ≥ 0.Then, every iteration leads to a KKT system of the following form, where slack increments Δs have already been eliminated: Here G, F, F r correspond to dynamic, global, and range constraints, respectively, and Φ ≥ 0, Ψ > 0 are diagonal barrier matrices.Elimination of Δv from the last equation, F r Δy + Ψ −1 Δv = −ψ, then yields a system of the form where If we solve the NLP by an active-set based SQP method, every active-set sub-iteration also produces a KKT system of the form (7), except that the barrier matrices vanish (yielding H = H ) and that F may contain additional rows from active inequality constraints.It is an essential feature of all three tree-sparse NLP forms that the additional terms in H (IPM) or in F (SQP) preserve the original block structure.
Next we consider the separable or partially separable Lagrangians of the three treesparse NLP forms and the resulting specializations of the KKT system (6), which are all straightforward extensions of the material in Steinbach (2002).The specific forms of the Lagrangians are needed for the Quasi-Newton updates in Sect. 3 whereas the KKT systems are needed for the inertia correction in Sect. 4.

Implicit tree-sparse KKT system
The Lagrangian of the implicit tree-sparse NLP (with where L j (y j , λ j , λ S( j) , v j , μ, η j ) = φ j (y j ) + λ T j g j (y j ) − k∈S( j) With barrier parameter β and e = (1, . . ., 1) T in suitable dimension, the resulting barrier matrices and vectors are Let ζ generically denote the dual variables appearing in the Lagrangian L or its components L j .Then the relevant partial derivatives of the Lagrangian read G j = ∇ y i h j (y i ), P j = ∇ y j g j (y j ), F j = ∇ y j e j (y j ), F r j = ∇ y j r j (y j ), and we further abbreviate This yields the following KKT system (6) in node-wise representation, where we omit Δ's on all increments and write y j , λ j , v j , μ for simplicity: The construction of H , G, F, F r and Φ, Ψ from the node contributions proceeds exactly as in Steinbach (2002), also for the following explicit tree-sparse NLP variants.

Outgoing tree-sparse KKT system
The Lagrangian of the outgoing tree-sparse NLP (with ξ j = (ξ − j , ξ + j ) and h 0 again constant) reads L(x, u, λ, v, μ, η, ξ) Here, the barrier matrices and vectors are and the partial derivatives of the Lagrangian read F r j = ∇ x j r j (x j , u j ), D r j = ∇ u j r j (x j , u j ), (11e) Thus, with we obtain the KKT system (6) in node-wise representation, again with Δ's omitted:

Incoming tree-sparse KKT system
For the incoming tree-sparse NLP the Lagrangian finally reads and Here, the barrier matrices and vectors are and the partial derivatives of the Lagrangian read Again, with we obtain the KKT system (6) in node-wise representation with Δ's omitted: 3 Tree-sparse Quasi-Newton updates Standard interior-point methods are second-order methods, i.e., they use local information given by the Hessian of the Lagrangian L. In many applications, the evaluation of second-order derivatives is prohibitively expensive and one is interested in using Quasi-Newton updates that approximate the Hessian.
For the following let an initial point y (0) and an initial approximation B (0) of the Hessian ∇ 2 yy L(y, ζ ) be given, where ζ again denotes all relevant dual variables.Then, omitting the current iteration index k and denoting iteration k + 1 by superscript "+", three standard approaches of updating the Hessian approximations are the symmetric rank-two BFGS update rule the symmetric rank-one (SR1) update formula and the somewhat less popular PSB rule where the iteration data is given by An overview of Quasi-Newton methods for unconstrained and constrained optimization can be found, e.g., in Dennis and Moré (1977), Dennis and Schnabel (1996), Fletcher (2013), Nocedal and Wright (2006) and the references therein.General-purpose implementations of algorithms using Quasi-Newton approaches for constrained nonlinear optimization include, e.g., KNITRO (Byrd et al. 2006(Byrd et al. , 2000) ) and SNOPT (Gill et al. 2002).

Quasi-Newton methods for tree-sparse problems
Standard update formulas like ( 16)-( 18) produce dense Hessian approximations that destroy the specific block structure of the tree-sparse NLPs.Quasi-Newton methods in large-scale optimization should generally not alter the sparsity pattern too much, i.e., the Hessian update strategy should keep additional fill-in at a minimum level.Such sparse Quasi-Newton approaches for unconstrained and constrained optimization are considered repeatedly in the literature; see, e.g., Fletcher (1995), Gill et al. (1984), Liu and Nocedal (1989), Lucia (1983), Powell and Toint (1981), Toint (1981).
In fact, the tree-sparse NLPs are designed such that their block structure is entirely preserved by suitable block-sparse Hessian approximations.The key idea is to apply the update formulas node-wise rather than globally, which goes back to multiple shooting SQP methods for deterministic optimal control problems (Bock and Plitt 1985).As an additional benefit the updates have much higher rank.
Before we explicitly discuss tree-sparse Quasi-Newton updates for specific control forms we briefly review the concept of (partially) separable functions.A function ϕ : R n → R is called partially separable if it can be written as where each function ϕ i depends only on a subset of the variables that is indexed by Quasi-Newton methods for arbitrary partially separable functions have first been developed in Griewank and Toint (1982a, b).

Tree-sparse Hessian updates for outgoing control problems
For tree-sparse NLPs in outgoing control form the Lagrangian (10), is completely separable with respect to the primal node variables, y j = (x j , u j ), j ∈ V .Thus, its Hessian is block-diagonal, and each block ∇ 2 y j y j L j (y j , ζ ) can be approximated individually by some B j obtained with the same update formula based on Here, B j approximates the diagonal block (without barrier terms) where the subblocks H j , J j , and K j are given in ( 11).The overall Hessian approximation (again without barrier terms) is then given by Finally note that, since the Lagrangian (8) of the implicit NLP form is also completely separable, it admits analogous individual updates of the diagonal blocks ∇ 2 y j y j L j (y j , ζ ) = H j .No communication is needed in the distributed code.

Tree-sparse Hessian updates for incoming control problems
The Lagrangian of tree-sparse NLPs in incoming control form, is composed of two types of node functions, L i j and L j as given in (13).In contrast to the outgoing control case, this Lagrangian is only partially separable with respect to the primal node variables y j = (u j , x j ).Thus, by setting y i j = (x i , u j ), the overall Hessian ∇ 2 yy L will be approximated by the sum of updates B i j ≈ ∇ 2 y i j y i j L i j and B j ≈ ∇ 2 x j x j L j , j ∈ V , whose subblocks overlap in part as detailed below.We obtain the node-wise approximations where H i j , J j , and K j are defined in ( 14).In a distributed code, B i j has to be sent to the process that holds B i .The approximations of the Hessians of L j are given by The approximations B i j and B j are updated by applying, for instance, one of the update rules ( 16)-( 18) where the respective differences of the iterates are the corresponding differences of the gradients of the Lagrangian are and, finally, r i j = g i j − B i j s i j and r j = g j − B j s j .
To illustrate the resulting block structure, we consider a simple tree with root 0 and successors S(0) = {1, 2}.In this case, the Hessian of the Lagrangian (13) has the form with variables (u 0 , x 0 , u 1 , x 1 , u 2 , x 2 ) and Hence, the subblocks K j and H j are placed on the diagonal corresponding to the respective variables u j and x j .Moreover, each subblock H i j is added to Hi , and J j and J T j are placed on the anti-diagonals corresponding to the variable pair (x i , u j ).

Tree-sparse inertia correction
Recall that the search direction in an interior-point method for NLP (1) is computed from system (7).Hence the reduced KKT matrix has to be invertible.In addition, to guarantee a descent direction, the Hessian H has to be positive definite on the nullspace of the constraints matrix.These conditions are satisfied if and only if the reduced KKT matrix has n positive and m negative eigenvalues (Nocedal and Wright 2006), i.e., inertia(Ω) = (n, m, 0).Therefore, when necessary, the inertia condition is enforced by replacing Ω with as follows.If C (and hence Ω) is rank-deficient, the associated zero eigenvalues are shifted into the negative region by choosing a fixed small value δ r > 0 (regularization).Then, if Ω(0, δ r ) has less than n positive eigenvalues, δ c is increased repeatedly until Ω(δ c , δ r ) has the desired inertia (convexification).For details see Schmidt ( 2013), Vanderbei and Shanno (1997), Wächter and Biegler (2006).
For the inertia correction of the tree-sparse NLPs, additional structural properties of the associated factorizations of Ω come into play.First, in the implicit tree-sparse KKT system the general regularization −δ r I in (27) splits into independent regularizations for each dynamics equation (9b) and for the global constraints (9d), which are applied to certain symmetric Schur complement blocks Y j ≥ 0 and X ∅ ≥ 0 (see Steinbach 2001) by factorizing Y j + δ r I > 0 (or X ∅ + δ r I > 0) unless Y j > 0 (or X ∅ > 0).Thus, any rank deficiencies in C are handled locally, and the local refactorizations do not require any communication in the distributed code.
On the other hand, the general regularization −δ r I would destroy the block structure of the tree-sparse KKT systems in outgoing or incoming control form since in both cases the factorization makes use of the fact that the matrix blocks to the right of G and below G T are zero.These factorizations can only handle a regularization of the global constraints, not the dynamics, While at first thought this might appear as a drawback, it is actually an advantage: with explicit dynamics (3b) or (4b), the term −x j creates identity blocks −I along the diagonal of G. Hence G has always full rank by construction, and a regularization needs only be considered for F. A closer look at the tree-sparse KKT solution algorithms reveals that full rank of F is equivalent to X ∅ > 0 (where X ∅ can be shown to be the Schur complement of the projection of F on N (G), with X ∅ ≥ 0).This is checked in the very last block operation: the Cholesky factorization of X ∅ , cf. step 18 of Table 1 for the incoming control case.Thus, if a regularization is required, only X ∅ + δ r I > 0 needs to be refactorized while the by far more expensive initial part of the factorization can be retained.Again there is no extra communication.Second, to preserve sparsity, each tree-sparse factorization needs a stricter condition than the inertia condition above: while rank(C) = m is required as before, H must be positive definite on the null space N (G) rather than on N (C) = N (G) ∩ N (F).This is because all three factorizations form the Schur complement X ∅ of the projection of F on N (G).The stricter condition can easily be incorporated into the convexification heuristic since positive definiteness of H + δ c I on N (G) is checked during the factorization.Of course, in problems without global constraints we have F ∈ R 0×m and there is no difference.
The structural details of the tree-sparse convexification are more involved than for the regularization.Positive definiteness of H + δ c I on N (G) is equivalent to positive In step 18, X ∅ , e ∅ denote X V , e V after finishing ("removing") every node j ∈ V in step 17.See Steinbach (2002) for details definiteness of properly modified blocks K j in every node, see steps 3 and 7 in Table 1.Thus, if this does not hold, the Cholesky factorization of some modified block K j will fail (step 8), and Ω * (δ c , δ r ) has to be refactorized from scratch with an increased value of δ c .In order to avoid refactorizations of the entire matrix, different convexification parameters δ cj can be used in each node j to modify the local Hessian blocks, This way, only local refactorizations are needed.However, in contrast to the uniform shift of all eigenvalues of H by δ c , individual local shifts δ cj may produce a highly inhomogeneous change of the subproblem's geometry.As a remedy, one may also consider combined global and local convexifications, as suggested in Hübner (2016).This allows for a wide range of possible heuristics that we do not want to explore here.
Let us finally discuss the inertia correction with local equality constraints.In all tree-sparse factorizations, these constraints are eliminated by local projections at the very beginning of the algorithm, under the assumption that the Jacobians of e l i j , e l j in (5) have full rank.This produces a projected KKT system without local constraints which has precisely the form ( 9), ( 12), or (15).If any of the local constraints may require a regularization, this procedure has to be modified as follows.All potentially rank-deficient local constraints are modeled as global constraints.However, it turns out that each of them augments the Schur complement X ∅ ≥ 0 mentioned above with an independent diagonal block X l j ≥ 0. Thus, similar to the implicit case, the general regularization splits into independent local regularizations that are applied immediately when reaching the respective block by factorizing X l j +δ r I > 0 (or X ∅ + δ r I > 0) unless X l j > 0 (or X ∅ > 0).Local regularizations are again automatically independent, require no communication, and in contrast to the convexifications a refactorization of Ω * (δ c , δ r ) from scratch is never needed.

Case studies
The main goal of this section is to show that the tree-sparse Hessian updates and the tree-sparse inertia correction work efficiently: they provide natural extensions of the existing tree-sparse algorithms to deal with nonconvexity and with unavailable or expensive second-order derivatives.
In Sect.5.1 and 5.2 we consider robust moving horizon control problems for a double integrator and a bioreactor, respectively.These examples are both nonconvex.The second problem does not provide explicit evaluations of the Hessian of the Lagrangian, so we use a Quasi-Newton approach based on the tree-sparse Hessian updates of Sect. 3.All optimization problems are solved using the interior-point code Clean::IPM (Schmidt 2013) with a tree-sparse KKT solver.The KKT solver incorporates the proposed inertia correction heuristic of Sect. 4 to address nonconvexity.
Moving horizon controllers (MHC) compensate random disturbances of a process by measuring or estimating the current disturbance in regular intervals and solving a dynamic optimization problem over a certain prediction horizon T to determine the current corrective action.Robust MHC incorporate an explicit stochastic model of future disturbances up to some stochastic horizon T s ≤ T whereas standard MHC simply ignore future disturbances (T s = 0).

Nonlinear double integrator
In this case study, we consider a moving horizon controller (MHC) for stabilizing a perturbed nonlinear double integrator.The dynamics is given by the discrete time model proposed in Lazar et al. (2008), Herein, the state (x 1 , x 2 ) is driven by the control u ∈ [−2, 2], and the first state x 1 is perturbed by some dynamic disturbance d with nominal value d nom = 0.The task is to keep the system close to its reference point, (x * 1 , x * 2 ) = (0, 0).In the absence of disturbances, the reference state (x * 1 , x * 2 ) is a fixed point of the dynamics (30).As in Lucia and Engell (2012), the disturbance may take one of three values, d(t) ∈ {−0.05, 0, 0.05}, with respective probabilities 0.2, 0.4, and 0.4.Again as in Lucia and Engell (2012), the objective is with Q = I and R = 0.15, a standard quadratic tracking functional with control costs.At each sampling time, the double integrator receives a control signal u(t) from the MHC and the random disturbance d(t) is observed.The resulting new state (x 1 (t + 1), x 2 (t + 1)) is then sent back to the controller; see Fig. 1.

Tree-sparse formulation
The moving horizon optimization problem that is solved at each sampling time to determine the control signal u includes the dynamic model ( 30) and the cost (31).We formulate it as a tree-sparse NLP in outgoing control form with probability-weighted node objectives and simple control bounds Fig. 2 Progress of average states of perturbed double integrator for three initial states with deterministic MHC (T s = 0, T = 3) The initial state x0 is given as h 0 = x0 .A possible scenario tree is shown in Fig. 1.

Control towards reference point
In a first test, we evaluate the performance of the deterministic controller (T s = 0) for bringing the system from three given initial states close to the reference point (x * 1 , x * 2 ).The MHC uses the prediction horizon T = 3 and analytic second-order derivatives.We generate random disturbances at t = 1, . . ., 10 and apply the MHC for 10 time steps.For each initial value, this is repeated 50 times with different series of random disturbances.Figure 2 shows the mean values of the states over the 10 time steps.For each of the three initial states, the mean values come very close to (x * 1 , x * 2 ) within the first 5 time steps, which confirms the proper operation of the controller.

Hold the reference point with minimal costs
In the second test series, the performance of the robust controller is tested for keeping the perturbed system close to the initial state (x * 1 , x * 2 ) over 20 time steps.Here we consider prediction horizons T = 3 and T = 10 with respective stochastic horizons T s ∈ {0, 1, 2, 3} and T s ∈ {0, 1, 2, 3, 5, 7}.Again we use analytic second-order derivatives, and each test is run with the same set of 50 series of random disturbances at t = 1, . . ., 20. Figure 3 illustrates the resulting averages of accumulated costs, states, and control over 20 steps.We observe that the 10 different controllers determined by (T , T s ) show visible differences in the mean values of u and x 2 only during the first five time steps (with an average control of zero at t = 5).Moreover, the length of the prediction horizon has no visible influence: the plots for T = 3 and T = 10 combined with T s ∈ {0, 1, 2, 3} look exactly identical.A close examination of the data reveals that differences after t = 5 and between T = 3 and T = 10 are indeed very small with values in the order of 10 −4 to 10 −3 .
The average accumulated costs (top of Fig. 3) demonstrate that the performance of the controller improves when uncertainties are included.Adding uncertainties in the first time step, i.e., increasing T s from 0 to 1, reduces the resulting costs at the final step t = 20 by approximately 9%.Increasing the stochastic horizon to T s = 2 leads to a further cost reduction of approximately 3%.Larger stochastic horizons (T s > 2) have no significant effect on the costs.The cost reductions are obtained by balancing the disturbances in the first 5 time steps of the process.The control signal is adjusted to reduce the deviation of x 1 from the reference state x * 1 = 0.This reduces the costs due to deviations of x 1 and increases the costs due to deviations of x 2 and due to nonzero control signals.In other words, within the critical first 5 time steps, the incurred costs are "moved" from state x 1 to state x 2 and to the control u.This strategy pays off in the further process where the costs due to u and x 2 are identical for all considered values of T s while the costs due to x 1 depend on its value at time step 5.

Exact Hessians versus approximations
Here we solve the robust control problem ( 32)-( 34) with a fixed prediction horizon T = 12 and increasing stochastic horizons T s ∈ {1, . . ., 12}.We measure the total time for IPM solution and for evaluating the NLP data with either analytic Hessians or SR1 updates or PSB updates.Each test is run 50 times for different initial states, and the averages of IPM solution time, NLP evaluation time, and iteration counts are measured, see Fig. 4. The problem sizes are given in Table 2 together with total runtimes of the IPM with analytic Hessian evaluations.All computations are carried out on a single core of a workstation with 48 GiB of RAM and 12 X5675 cores running at 3.07 GHz.
Let us first consider the case of analytic Hessian evaluations.The average number of iterations is almost independent of T s , which indicates a certain "scalability" of the stochastic model.The NLP evaluation time is lower than with both Quasi-Newton updates, thus analytic Hessians are cheaper than approximations in this problem.Finally we note that the constraints matrices cannot become rank-deficient and thus the regularization remains inactive for all 12 values of T s and all 50 initial values.
With SR1 updates we observe the lowest iteration counts of all three Hessian versions on small problems (T s < 8) but the highest iteration counts on large problems (T s > 9).The convexification is again inactive in all runs, which indicates that the tree-sparse SR1 updates provide good Hessian approximations.
With PSB updates, the iteration count varies with T s but remains significantly smaller than with analytic Hessians.This yields the smallest total runtimes on large problems (T s > 9) although the NLP evaluation time per iteration is significantly larger than for analytic Hessians.Surprisingly, most runs require a convexification to solve the optimization problems.

Nonlinear bioreactor
This case study is concerned with a robust nonlinear moving horizon controller that keeps a bioreactor in a steady state of production.The problem is proposed as a benchmark in Ungar (1990) and has been studied, e.g., in Lucia and Engell (2012), Lucia et al. (2012).The plant consists of a continuous flow stirred tank reactor containing a mixture of water and cells.The latter consume nutrients and produce (desired and undesired) products as well as more cells.The volume of the mixture is constant and its composition is adjusted by a water stream that feeds nutrients into the tank at the inlet and that contains nutrients and cells at the outlet.The dynamic model of the bioreactor is given by the ODE system where x 1 and x 2 are dimensionless amounts of cell mass and nutrient, respectively, and v is the flow rate of the water stream.The respective physical bounds are The ODE system (35) describes the rates of change in the amounts of cells x 1 and nutrients x 2 , respectively, that result from the respective amounts −x 1 v and −x 2 v leaving the tank and from the metabolism of the cells.The cell growth is represented by x 1 (1 − x 2 )e x 2 /γ , where γ is the uncertain nutrient consumption parameter with nominal value γ nom = 0.48.The rate of cell growth β depends very mildly on the composition of the mixture in the tank.We regard it as constant, β nom = 0.02.When feeding nutrients to the bioreactor with a constant flow rate v, the system has a Hopf bifurcation at a certain flow rate v H that depends on the values of γ and β.For v < v H , System (35) stabilizes at a unique fixed point (x * 1 , x * 2 ) whereas it becomes unstable for v ≥ v H .For the nominal parameter values, the Hopf bifurcation occurs at the flow rate v H = 0.829.This value decreases with an increasing value of γ or a decreasing value of β as shown in Fig. 5.
The desired steady state of production, (x * 1 , x * 2 ) ≈ (0.1236477, 0.8760318), is close to the Hopf bifurcation and is obtained for the constant flow rate v * = 0.769 with nominal parameter values.The parameter γ is assumed to be normally distributed with small variance, γ ∼ N (γ nom , 0.005).As in Lucia and Engell (2012), standard quadratic costs are applied that penalize deviations from the reference state x * 1 with the factor 200 and changes of the flow rate with the factor 75.

Tree-sparse formulation
The optimization problem for the bioreactor is formulated as a tree-sparse NLP with incoming control (4).State variables x j ∈ R 3 consist of the two states x 1 , x 2 , and the flow rate v of the ODE system (35).The control u j ∈ R models the change of the piecewise constant flow rate v at time t(i), yielding dynamics x j = g j (x i , u j ) = (w j (x i , u j ), x i,3 + u j ) where Here the flow rate is modeled as a state variable since its difference u j = x j,3 − x i,3 enters into the objective, with R = 75, Q = Diag(200, 0, 0), and x * = (x * 1 , x * 2 , 0).Finally, physical bounds (36) are incorporated as simple state bounds, x − j = (0, 0, 0) and It turns out that with increasing stochastic horizon T s the robust control problem becomes extremely hard to solve due to the highly nonlinear dynamics.Therefore we use as scenario tree a simple fan (T s = 1) with n s scenarios, see Fig. 6.

Keep a steady state of production
In the following, we consider 13 instances of the benchmark problem differing in the prediction horizon T or in the number of scenarios n s .The reactor is started in the reference state and run for 40 s with a sampling interval of 100 ms.At each sampling time, the plant receives a new control signal from the MHC and the random value of γ is sampled, see Fig. 6.Thus, each test run requires solving 400 optimization problems.Each of those tests is run 50 times with different series of random disturbances.We compute variances of the cell mass x 1 and averages of accumulated costs, of x 1 , and of the flow rate v.
As in Sect.5.1, the tree-sparse NLPs are solved using Clean::IPM where solutions of the IVP (37) are computed using a semi-implicit extrapolation method (Bader and Deuflhard 1983).The Hessians are approximated by tree-sparse SR1 updates, which turned out to be the most reliable choice here.All computations in this section have been executed on a single core of a workstation with 16 GiB of RAM and an Intel(R) Core i7-3770 quad-core processor running at 3.40 GHz.
In some test runs, the IPM does not converge for all 400 optimization problems, which demonstrates the difficulty of these problems.To keep the results comparable, we exclude all test runs that do not succeed in each of the 13 benchmark instances.In the end we are left with solutions for 28 of the 50 disturbance series.
The results are given in Figs.7 and 8.The plots show averages of accumulated costs (upper left), of the cell mass x 1 (upper right), the variance of x 1 (lower left), and average flow rates (lower right) vs. time.Two situations are shown in the figures.First, in Fig. 7, the bioreactor is regulated using deterministic optimization problems (n s = 1) and an increasing prediction horizon T .Second, in Fig. 8, the MHC uses stochastic optimization problems with fixed prediction horizon (T = 5) and a varying number of scenarios n s .In both situations, increasing the free tree parameter (T or n s ) improves the performance of the controller.The average accumulated costs decrease monotonously, which is achieved by damping the effect of the perturbation on the cell mass x 1 .After experiencing disturbances, x 1 oscillates about the reference state x * 1 .The amplitude of this oscillation is reduced by adjusting the flow rate v, which reduces the costs caused by the deviations of x 1 .
The results in Fig. 7 show that a certain minimal prediction horizon is required to produce a significant effect.Using the value T = 5 reduces the costs by less than 5%  in comparison to T = 1 whereas the value T = 10 yields a reduction of more than 60%.However, increasing T further yields no significant improvement.The value T = 20, for instance, produces a reduction of only 65%.The results in Fig. 8 show that the robust control approach with its explicit model of uncertainties improves the controller performance.The total costs are reduced by 7% with n s = 3, and by 70% with n s = 81.Of course, each additional scenario causes higher computational effort while the relative reduction of costs per scenario decreases.For instance, we have 3.5% reduction of costs per scenario for n s = 3 and 0.88% per scenario for n s = 81.
The results of Figs. 7 and 8 show that the goal of reducing the costs is achieved in both situations.However, increasing the number of scenarios n s is computationally more expensive than enlarging the prediction horizon T , as can be seen in Table 3.The optimization problems of case 1 (T = 100 and n s = 1) are significantly smaller and solved faster than those of case 2 (T = 5 and n s = 51), although the two cases produce almost identical costs as shown in Fig. 9. Finally, case 3 in Table 3 represents a fair compromise between the length of the prediction horizon (T = 10) and the number of scenarios (n s = 3).In Fig. 9 we see that this compromise produces the smallest costs and, moreover, has the smallest optimization problems among the ones listed in Table 3.It is solved 3 times faster than case 1 and 9 times faster than case 2. Thus we conclude that a suitably designed robust controller can be significantly more efficient than a deterministic controller.

Conclusion
We have seen that the concept of tree-sparse optimization problems (Steinbach 2002) extends directly from convex problems with polyhedral constraints to the fully non-linear and nonconvex NLP case.The required structural assumptions are satisfied whenever the nonlinear constraint functions possess mild separability properties.Moreover, interior-point methods with a tree-sparse KKT solver are directly applicable to the general case by introducing straightforward structure-specific specializations to a standard inertia correction procedure.Finally, the Lagrangian is separable or partially separable in all three tree-sparse NLP variants so that the Hessian can be approximated by structure-preserving Quasi-Newton update formulas when second-order derivatives are unavailable or too expensive.The effectiveness of these techniques has been demonstrated on two robust model predictive control problems with ODE dynamics.Altogether, our problem formulation and our solution approach exploit the specific structure of tree-sparse problems in a natural way.
Although we have not demonstrated this here, it is easily seen that the distributed implementation presented in Hübner et al. (2017) extends as well to the NLP case considered in this paper: the function and derivative evaluations, possibly with Quasi-Newton updates, and the tree-sparse inertia correction parallelize similar or easier than the KKT solver.Our techniques are also directly applicable in an active set-based SQP framework such as Rose ( 2018): here we can again use the tree-sparse KKT solver, except that the implementation requires some overhead for the changing active set.

Fig. 1
Fig. 1 Moving horizon control of the double integrator (left) and scenario tree with prediction horizon T = 3 and stochastic horizon T s = 2 (right)

Fig. 4
Fig. 4 Double integrator: average IPM performance vs.stochastic horizon T s for different algorithm configurations

Fig. 5
Fig. 5 Hopf bifurcations of the bioreactor for varying parameters γ and β

Fig. 6
Fig. 6 Process control of the bioreactor (left) and scenario tree with prediction horizon T = 3 and n s = 5 scenarios (right)

Fig. 7 Fig. 8
Fig. 7Performance of the bioreactor for deterministic problems (n s = 1) with increasing prediction horizons (T )

Table 1
Operations of KKT solution algorithm for incoming control

Table 2
Double integrator: problem sizes and average IPM solution times with analytic Hessian evaluations

Table 3
Problem sizes, IPM solution times (s) and NLP evaluation times (s) of the problems in Fig.9