A Numerical Method for Nonconvex Multi-objective Optimal Control ...
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems C. Yal¸cın Kaya∗
Helmut Maurer†
September 3, 2013
Abstract A numerical method is proposed for constructing an approximation of the Pareto front of nonconvex multi-objective optimal control problems. First, a suitable scalarization technique is employed for the multi-objective optimal control problem. Then by using a grid of scalarization parameter values, i.e., a grid of weights, a sequence of single-objective optimal control problems are solved to obtain points which are spread over the Pareto front. The technique is illustrated on problems involving tumor anti-angiogenesis and a fed-batch bioreactor, which exhibit bang–bang, singular and boundary types of optimal control. We illustrate that the Bolza form, the traditional scalarization in optimal control, fails to represent all the compromise, i.e., Pareto optimal, solutions.
Key words: Multi-objective optimization, Optimal control, Pareto front, Numerical methods, Tumor anti-angiogenesis, Fed-batch bioreactor.
1
Introduction
In optimal control, the task is to minimize a functional, referred to as the cost, typically stated in the so-called Bolza form, tf f0 (x(t), u(t), t) dt , (1) w ψ(x(tf ), tf ) + (1 − w) 0
subject to differential and algebraic constraints, where x(t) and u(t) are the state and control variables, respectively, and w is some fixed weight selected in the interval (0, 1). The terminal time tf might be fixed or free. The terminal cost, ψ(x(tf ), tf ), is usually in conflict with the integral cost in the second term; that is why a compromise solution has to be found depending on the value of the weight w. In the multi-objective (or vector) optimization terminology, the Bolza form in (1) is the so-called weighted-sum scalarization of the multi-objective problem with the two-objective functional vector, tf f0 (x(t), u(t), t) dt , (2) ψ(x(tf ), tf ), 0
to minimize (in the Pareto sense), subject to the same constraints as those of the Bolza optimal control problem with (1). ∗
School of Information Technology and Mathematical Sciences, University of South Australia, Mawson Lakes, S.A. 5095, Australia. E-mail: [email protected] . † Institut f¨ ur Numerische und Angewandte Mathematik, Westf¨ alische Wilhelms-Universit¨ at M¨ unster, M¨ unster, Germany. E-mail: [email protected] .
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
2
Broadly speaking, the Pareto minimization of (2) is the process of finding a compromise solution, referred to as a Pareto minimum, where the value of the terminal cost or the value of the integral cost cannot be improved (i.e., reduced) further, without making the value of the other cost worse (i.e., higher). The set of all such compromise solutions form the Pareto set, or the Pareto front. While one would expect to obtain all Pareto minima by solving the single-objective Bolza problem for all w ∈ [0, 1], it turns out that this may not be possible. In other words, it does not seem to be a good idea to use the Bolza form (1) of the objective functional for nonconvex optimal control problems, as solutions found over the range of weights may not yield all of the compromise solutions. In this paper, we propose to use a scalarization different from the Bolza, or the weightedsum, form for the terminal and integral costs: we combine these costs by means of a max operator and weights, which is referred to as the weighted Tchebychev scalarization in the literature [33]. This kind of scalarization is certainly not new for general multi-objective optimization problems, see, for example, [14, 33], in the finite dimensional setting; however, to the authors’ knowledge it has not been used for optimal control problems. The weighted Tchebychev scalarization is surjective from the space of weights to the Pareto set (or Pareto front). In other words, over the range of weights, the solution of the scalarized problem yields all of the Pareto minima. Conditions of optimality and related theoretical aspects of multi-objective optimal control problems have been studied in the literature rather extensively – see for example [10, 11, 17, 24]. Moreover, effective numerical methods have been developed and illustrated for convex problems – see [3, 6] and the references therein. However, there does not seem to be many numerical studies dealing with nonconvex multi-objective optimal control problems; among existing studies we note [28, 29, 30]. In these references, first the multi-objective control problem is discretized, and then the normal boundary intersection (NBI) and normalized normal constraint (NNC) techniques are employed as a scalarization of the resulting finitedimensional problems. This process can be summarized as discretize-scalarize-then-optimize. It is worth mentioning the work in [34], where also a discretize-scalarize-then-optimize approach is adopted as in [28, 29, 30], however with the main difference that the discretization is carried out in the context of mechanics involving unconstrained dynamical equations. Optimal control problems are infinite-dimensional optimization problems, and as such the scalarization should be defined and justified for the original infinite-dimensional problem. In this paper, we choose the weighted Tchebychev scalarization for the original problem, and then consider the scalarized problem’s discretization. Therefore, our approach can be summarized as scalarize-discretize-then-optimize. In any scalarization scheme based on weights, it is hard to know the meaning of the weights used for the objective functions. It is often not even possible to know the range of values of the weights over which the solutions change qualitatively. An advantage of the weighted Tchebychev scalarization is that it is possible to determine a subinterval of the “meaningful” weights, outside of which the Pareto solutions are the same as those obtained by utilizing the weights at the boundary of the subinterval. This is a particularly useful feature (a feature which is not possessed by the weighted-sum scalarization, i.e., the Bolza form), as then the computational resources can be utilized only for finding (different) Pareto solutions corresponding to the weights in the (meaningful) subinterval mentioned above. In the NBI and NNC techniques employed in [28, 29, 30], the meaningful interval of weights would be the whole interval [0, 1]. However, the weights that are used are not associated directly with the weights of the objective functionals we consider. In the current paper, we present an algorithm for generating an approximation of the
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
3
Pareto front over the meaningful subinterval of weights. We illustrate the working of the algorithm by means of two computationally challenging problems, one involving tumor antiangiogenesis from [25] and the other a fed-batch bioreactor from [29]. For these problems, the Pareto fronts are constructed and solutions corresponding to different weights in the Pareto front are displayed. The first problem results in solutions with bang–bang and singular type of controls. The second problem involves a richer structure: the solutions, depending on the values of the weights, contain bang, singular and boundary arcs, the latter being due to the presence of (active) pure state constraints. We numerically verify that the solutions we find indeed satisfy the necessary conditions of optimality. Note that the multi-objective fed-batch bioreactor problem is also solved in [29], but over a much coarser time-grid, and the solutions found are not checked to see whether or not they satisfy the necessary conditions of optimality. We solve the multi-objective tumor anti-angiogenesis problem for the first time, in this paper. The paper is organized as follows. In Section 2, we describe the multi-objective optimal control problem, and provide definitions pertaining to scalarization. We also cite a result, which can readily be extended to the optimal control setting we have, ensuring that the Tchebychev scalarization we employ will be successful. In Section 3, we provide the algorithm. In Section 4, we implement the algorithm on the two practical optimal control problems and discuss the subsequent numerical results. In Section 5, we make concluding remarks and identify some future work.
2
Problem Statement and Preliminaries
Consider the multi-objective optimal control problem min (ϕ1 (x(tf ), tf ), . . . , ϕr (x(tf ), tf )) ˙ = f (x(t), u(t), t) , for a.e. t ∈ [0, tf ] , subject to x(t) (OCPr ) θ(x(0), x(tf ), tf ) = 0 , C(x(t), u(t), t) ≤ 0 , for a.e. t ∈ [0, tf ] , S(x(t), t) ≤ 0 , for all t ∈ [0, tf ] , where the state variable x ∈ W 1,∞ (0, tf ; IRn ), x˙ := dx/dt, and the control variable u ∈ L∞ (0, tf ; IRm ). The functions ϕi : IRn × IR+ → IR, f : IRn × IRm → IRn , θ : IRn → IRp1 , C : IRn × IRm → IRp2 , and S : IRn → IRp3 , are continuous in their arguments. In this problem, tf can be fixed or free. Here, L∞ (0, tf : IRn ) corresponds to the space of essentially bounded, measurable functions equipped with the essential supremum norm. Furthermore, W 1,∞ (0, tf ; Rn ) is the Sobolev space consisting of functions x : [0, tf ] → Rn whose first derivatives lie in L∞ . Assume that ϕi (x(tf ), tf ) ≥ 0, for all i = 1, . . . , r. Note that this assumption can easily be met by adding a large enough positive number to each objective functional. Note that Problem (OCPr ) is in general a nonsmooth problem, because it does not require differentiability of the objective functionals or the constraints. Moreover, although we have stated Problem (OCPr ) in very broad terms, it can further be generalized, for example by adding multi-point constraints, partial differential equations, etc. In other words, although Problem (OCPr ) is already in a more general form than what one usually encounters in applications, it can be further made look more general. , where tmax > 0 is some constant. Define the For technical convenience, let tf ≤ tmax f f
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
4
feasible set, X ⊂ W 1,∞ (0, tf ; IRn ) × L∞ (0, tf ; IRn ) × IR+ , such that ˙ = f (x(t), u(t), t) , for a.e. t ∈ [0, tf ]; θ(x(0), x(tf ), tf ) = 0 ; X := {(x, u, tf ) : x(t) C(x(t), u(t), t) ≤ 0 , for a.e. t ∈ [0, tf ]; S(x(t), t) ≤ 0, for all t ∈ [0, tf ]} . Define the vector of objective functionals, ϕ(x(tf ), tf ) := (ϕ1 (x(tf ), tf ), . . . , ϕr (x(tf ), tf )). The triplet (x∗ , u∗ , t∗f ) ∈ X is said to be a Pareto minimum if there exists no (x, u, tf ) ∈ X such that ϕ(x(tf ), tf ) = ϕ(x∗ (t∗f ), t∗f ) and ϕi (x(tf ), tf )) ≤ ϕi (x∗ (t∗f ), t∗f ) ,
for all i = 1, . . . , r .
On the other hand, (x∗ , u∗ , t∗f ) ∈ X is said to be a weak Pareto minimum if there exists no (x, u, tf ) ∈ X such that ϕi (x(tf ), tf )) < ϕi (x∗ (t∗f ), t∗f ) ,
for all i = 1, . . . , r .
The set of all vectors of objective functional values at the Pareto and weak Pareto minima is said to be the Pareto front (or efficient set) of Problem (OCPr ) in the r-dimensional objective value, or outcome, space. Note that the coordinates of a point in the Pareto front are simply ϕi (x∗ (t∗f ), t∗f ), i = 1, . . . , r. Obviously, when r = 2, the Pareto front is in general a curve; and when r = 3, the Pareto front is in general a surface.
2.1
Scalarization
For computing a solution of the nonconvex multi-objective optimal control problem (OCPr ), we will consider the following single-objective problem, i.e., scalarization. This scalarization is often considered for nonconvex multi-objective finite-dimensional optimization problems [33]. (Pw )
min
(x,u,tf )∈X
max{w1 ϕ1 (x(tf ), tf ), . . . , wr ϕr (x(tf ), tf )} ,
where wi , i = 1, . . . , r, are referred to as weights. The following theorem is a direct consequence of [20, Corollary 5.35]. It can also be concluded by generalizing the proofs given in a finite-dimensional setting in [33, Theorems 3.4.2 and 3.4.5] to the infinite dimensional setting in a straightforward manner. Theorem 1 The triplet (x∗ , u∗ , t∗f ) is a weak Pareto minimum of (OCPr ) if, and only if, (x∗ , u∗ , t∗f ) is a solution of (Pw ) for some w1 , . . . , wr > 0. An ideal cost ϕ∗i , i = 1, . . . , p, associated with Problem (Pw ) is the optimal value of the optimal control problem, min ϕi (x(tf ), tf ) . (Pi ) (x,u,tf )∈X
Let (x, u, tf ) be a minimizer of the single-objective problem (Pi ) above. Then we also define ϕj := ϕj (x(tf ), tf ), for j = i and j = 1, . . . , r. Problem (Pw ) is referred to as the weighted Tchebychev problem (or Tchebychev scalarization) because of the weighted Tchebychev norm, maxi |wi ϕi (x(tf ), tf )| = maxi wi ϕi (x(tf ), tf ), appearing in the objective. In the case when ϕ∗i is negative, one can simply add a large enough positive number to the ith objective, to make the objective positive. In general, it is useful to add a positive number to each objective in order to obtain an even spread of the Pareto points approximating the
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
5
Pareto front – see, for example, [14], for further discussion and geometric illustration. For this purpose, we define next the so-called utopian objective values. A utopian objective vector β ∗ associated with Problem (OCPr ) consists of the components given as βi∗ = ϕ∗i − εi where εi > 0 for all i = 1, . . . , p. Problem (Pw ) can be equivalently written as βi∗
min
(x,u,tf )∈X
max{w1 (ϕ1 (x(tf ), tf ) − β1∗ ), . . . , wr (ϕr (x(tf ), tf ) − βr∗ )} .
In the case when the objective functionals and the constraints in Problem (OCPr ) are differentiable in their arguments, Problem (Pw ), or its equivalent above, is still a non-smooth optimization problem, because of the max operator in the objective. However, it is well-known that Problem (Pw ) can be re-formulated in the following (smooth) form. min α α≥0 (x,u,t )∈X f subject to w1 (ϕ1 (x(tf ), tf ) − β1∗ ) ≤ α , (OCPw ) .. . wr (ϕr (x(tf ), tf ) − βr∗ ) ≤ α . Problem (OCPw ) is referred to as goal attainment method [33], as well as Pascoletti-Serafini scalarization [15]. We will solve Problem (OCPw ) in an algorithm we present in the next section, for the two examples we want to study. Note that the “popular” weighted-sum scalarization, given below, fails to generate the nonconvex parts of a Pareto front. (Pws )
min
(x,u,tf )∈X
r
wi ϕi (x(tf ), tf ) .
i=1
We will illustrate this in the fed-batch bioreactor problem in Section 4.2. It should be noted that, for r = 2, the above weighted-sum scalarization can be shown to be equivalent to the Bolza form in (1), after defining a new state variable, as in the case of the tumor antiangiogenesis problem in Section 4.1.
3
An Algorithm to Generate the Pareto Front
In this section, we provide an algorithm which generates an approximation of the Pareto front of Problem (OCPr ) with two objectives. Note that by choosing w1 = w, and w2 = 1 − w, where w ∈ [0, 1], one can simply consider the single weight w in the scalarized problem (OCPw ). Solution of Problem (OCPw ) using equal increments of the value of w in [0, 1] can generate a reasonable approximation of the Pareto front, if one chooses a utopia point (β1∗ , β2∗ ) far enough from the front (see [14]). With the Tchebychev scalarization, the weight w would usually take values over a (smaller) subinterval of [0, 1]. In other words, a restricted range of values of w, i.e., 0 < w0 < w < wf < 1 would usually suffice for the generation of the front. Figure 3 illustrates the geometry with which one can compute the subinterval end-points, w0 and wf . Note that the points (ϕ∗1 , ϕ2 ) and (ϕ1 , ϕ∗2 ) represent the boundary of the Pareto front, where the equations of the “rays” emanating from the utopia point (β1∗ , β2∗ ) and passing through the boundary points are given as in the figure. Also note that these equations are obtained by virtue of the weighted
6
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
ϕ2
wf (ϕ1 − β1∗ ) = (1 − wf ) (ϕ2 − β2∗ ) (ϕ∗1 , ϕ2 ) Pareto front w0 (ϕ1 − β1∗ ) = (1 − w0 ) (ϕ2 − β2∗ ) (ϕ1 , ϕ∗2 ) ϕ1
(β1∗ , β2∗ )
Figure 1: Determination of the boundary weights w0 and wf .
Tchebychev (or max) norm. By substituting the boundary values of the Pareto curve into the respective equations, and solving each equation for w0 and wf one simply gets w0 =
(ϕ∗2 − β2∗ ) (ϕ1 − β1∗ ) + (ϕ∗2 − β2∗ )
and
wf =
(ϕ2 − β2∗ ) . (ϕ∗1 − β1∗ ) + (ϕ2 − β2∗ )
(3)
Although the algorithm below is described for the bi-objective case (i.e., when r = 2), it can be generalized to problems with more than two objectives, by re-defining the boundary of the space of weights and generating a grid of weights on a segment of a hyperplane rather than a subinterval. However, care needs be taken in doing this: for example for r = 3, it is well known that the boundary of the Pareto front (which now typically is a closed curve rather than just two points) may very well lie outside the triangle formed by the three points, (ϕ∗1 , ϕ2 , ϕ3 ), (ϕ1 , ϕ∗2 , ϕ3 ) and (ϕ1 , ϕ2 , ϕ∗3 ). While this kind of generalization is quite important, it is outside the scope of the present paper. Algorithm 1 Step 0.0 (Initialization) Choose the utopia parameters, ε1 , ε2 > 0. Set the number of discrete (approximating) points, (N + 1), in the Pareto front. Set k := 1. Step 0.1 (Boundary of the front) 1
1
1
2
2
2
(a) Find (x1 , u1 , tf ) that solves Problem (P1 ). Note that ϕ∗1 = ϕ1 (x1 (tf ), tf ) and 1 1 ϕ2 = ϕ2 (x1 (tf ), tf ). Mark a boundary point in the Pareto front: ϕN := (ϕ∗1 , ϕ2 ) . (b) Find (x2 , u2 , tf ) that solves Problem (P2 ). Note that ϕ1 = ϕ1 (x2 (tf ), tf ) and 2 2 ϕ∗2 = ϕ2 (x2 (tf ), tf ). Mark a boundary point in the Pareto front: ϕ0 := (ϕ1 , ϕ∗2 ) .
Step 0.2 (Utopia point) Set β ∗ := (β1∗ , β2∗ ) with βi∗ := ϕ∗i − εi , i = 1, 2. Step 0.3 (Range of weights) Determine the range of (meaningful) weights, i.e., the subinterval [w0 , wf ], using (3). Set the increment, ∆w := (wf − w0 )/N .
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
7
Step k.1 (Current weights) Set w := w0 + k ∆w. Set w1 := w and w2 := 1 − w. Step k.2 (A Pareto minimum) Find (x∗ , u∗ , t∗f ) that solves Problem (OCPw ). Assign a point in the Pareto front: ϕk := (ϕ1 (x∗ (t∗f ), t∗f ), ϕ2 (x∗ (t∗f ), t∗f )). Step k.3 (Stopping criterion) If k = N then STOP. Otherwise, set k := k + 1 , and go to Step k.1. In Step k.1 of Algorithm 1, the weights are chosen equidistantly. One should note that for a uniform distribution of the points in the value space approximating the Pareto front, these weights can as well be chosen in an adaptive fashion as discussed by Eichfelder in [15].
4
Numerical Examples
In this section, we illustrate the working of the algorithm on two numerical example problems, one involving tumor anti-angiogenesis in Section 4.1 and the other a fed-batch bioreactor in Section 4.2. It is common practice to approximate the state and control variables over a partition of the time horizon in optimal control problems, a process referred to as discretization, so that one can solve a (large-scale) finite-dimensional optimization problem, instead of an infinitedimensional one, to obtain an approximate solution – see, for example, [4, 5, 12, 13, 18, 21, 22, 31]. This overall process is referred to as a direct method, or a discretize-then-optimize approach. Convergence theory for the discretize-then-optimize approach and its applications have been studied in the literature extensively; see the above-listed references, which in general concern Runge-Kutta (and, as a particular case, Euler) discretization of state- and controlconstrained optimal control problems. In [12], where Euler discretization is considered, for example, under the rather general assumptions that the control variable is continuous and second-order sufficient conditions hold, the discrete (i.e., approximate) solutions of the state, control and adjoint variables converge to the continuous-time solutions (the state and costate variables with respect to the W 1,∞ -norm and the control variable with respect to the L∞ norm). Higher-order Runge-Kutta discretization schemes, which require higher-order differentiability of the solution, in general provide higher-order convergence rates [18]. In the fed-batch bioreactor example, we discretize the problem over a sufficiently fine grid. The solution for the control variable is discontinuous, incorporating bang–bang, singular and boundary arcs. So far, there exists no formal proof of convergence that the discretized solution for this example converges to the infinite-dimensional solution in topology of the L1 space. Such a convergence result has been obtained by Alt et al. [1, 2] only in the case of a linear dynamics and a bang–bang solution. A convergence result for more general problems involving discontinuous solutions has not been established yet in the literature. Nevertheless, our discretization approach clearly shows the control structure of the fed-batch bioreactor example. In the tumor anti-angiogenesis example, the solution also involves bang–bang and singular controls. However, since we can derive a feedback expression for the candidate singular optimal control, we parameterize the problem with respect to the (unknown) times when the bang–bang and singular controls switch, and so we reduce the scalarized problem into what is called an Induced Problem, which is explained in some detail in Section 4.1. We still employ a discretization scheme for computing the bang-bang and singular arcs precisely enough.
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
8
In Steps 0.1 and k.2 of the algorithm for the example problems that we study, a discretized version of Problem (OCPw ) is solved by using Ipopt, version 3.2.4s, a popular optimization software based on an interior point method; see [36]. We use AMPL [16] as an optimization modelling language, which employs Ipopt as a solver. Note that, by Theorem 1, there is a bijection between the space of weights and the space of Pareto minima. So, as the grid of weights is refined, one would expect to find more points in the Pareto front and obtain a better approximation, ultimately converging to the true Pareto front. Recall that here we take the scalarize-discretize-then-optimize approach, meaning that we first scalarize the multi-objective problem in the infinite-dimensional space, and then discretize before we apply a finite-dimensional optimization method to find a discrete approximate solution of the scalarized problem. By the existing theory of discretization we outlined above, under certain assumptions, one can guarantee that (or at least check the necessary and/or sufficient optimality conditions a posteriori to observe that) the discrete approximate solution converges to a solution of the continuous-time scalarization of the original problem, yielding a Pareto minimum of the original problem. On the other hand, if one takes the existing discretize-scalarize-then-optimize approach, then scalarization is carried out for the discretization of the original multi-objective problem. To the authors’ knowledge, no theory is available for guaranteeing convergence of a solution of the discretized multi-objective control problem to a solution of the continuous-time multiobjective control problem.
4.1
Optimal protocols for tumor anti-angiogenesis
Let us cite from the Introduction in the article of Ledzewicz, Maurer and Sch¨attler [25]: “Tumor anti-angiogenesis is an indirect cancer treatment approach with the aim to limit the tumor’s growth, possibly even shrink the tumor, by depriving it of the vasculature it needs for a steady supply with nutrients.” In a seminal paper, Ledzewicz and Sch¨attler [27] formulated the anti-angiogenic cancer treatment as an optimal control problem and identified the optimal control structure as a bang–singular–bang control. Here, we consider the dynamical model for anti-angiogenic monotherapy in [25] which is slightly different from the one in [27]. The state and control variables of the problem are as follows. p : primary tumor volume [mm3 ], q : carrying capacity, or endothelial support [mm3 ], u : anti-angiogenic agent. We introduce two competing objective functionals to minimize : (i) The final tumor volume : p(tf ) , (ii) The final tumor volume, plus a factor of the total amount of anti-angiogenic agent administered : tf u(t) dt . p(tf ) + 140 0
The factor 140 in the second functional has been employed to guarantee the same optimal control structure for all weights of interest. Note that the terminal time tf is free. It may t at first seem that it is “more natural” to consider 0 f u(t) dt alone as the second objective functional; however, the minimization of this integral cost (for finding a boundary point of
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
9
the front) yields the trivial solution u(t) = 0 and tf = 0 (tf is free), which is not realistic. Therefore we consider a combination of p(tf ) and the integral cost with a large positive coefficient. Incorporating the dynamical equations and constraints as described following bi-objective optimal control problem: min (p(tf ), p(tf ) + 140 y(tf )) p subject to p˙ = −0.084 p ln q , (E1) q˙ = 5.85 q 2/3 − 0.00873 q 4/3 − 0.02 q − 0.15 q u , y˙ = u , y(tf ) ≤ 45 , 0 ≤ u(t) ≤ 15 .
in [25], we get the
p(0) = 8000 , q(0) = 10000 , y(0) = 0 ,
The scalarization of Problem (E1) in the form of problem (OCPw ) in Section 2 yields the following optimal control problem. min α subject to the same constraints as in (E1) and (E1w ) w (p(tf ) − β1∗ ) ≤ α , (1 − w) (p(tf ) + 140 y(tf ) − β2∗ ) ≤ α . Note that, in the above problem, ϕ2 ≥ ϕ1 . This has a geometric implication on the Pareto front and the weights that are associated in constructing it: the front simply lies above the line ϕ2 = ϕ1 in the ϕ1 -ϕ2 plane, and the meaningful interval of weights is, for β ∗ = (0, 0), a subinterval of [0.5, 1]. Algorithm 1, with N = 50 and β ∗ = (0, 0), results in the approximation of the Pareto front shown in Figure 2(a). The “meaningful” interval of weights is [0.628, 0.792], which is a subset of [0.5, 1] as expected. The solution for all w ∈ [0, 0.628] is the same, corresponding to one boundary point in the Pareto front, so is the solution for all w ∈ [0.792, 1], corresponding to the other boundary point. The switching function σu (t) = −0.15 λq (t) q(t) + λy (t) , where λq and λy are the adjoint (or costate) variables associated with the differential equations for the state variables q and y, determines the minimizing control according to if σu (t) < 0 15 , (4) if σu (t) > 0 u(t) = 0 , singular , if σu (t) = 0 ∀ t ∈ [t1 , t2 ] ⊂ [0, tf ] ¨u (t) = 0 ∀ t ∈ [t1 , t2 ], A singular control can be computed from the relations σu (t) = σ˙ u (t) = σ since any singular arc has order one. In [25, 26], the following feedback expression was obtained for a singular control:
5.85 − 0.00873 q 2/3 1 5.85 + 0.00873 q 2/3 + 3 (0.084) − 0.02 . using (p, q) = 0.02 q 1/3 5.85 − 0.00873 q 2/3 Moreover, all singular trajectories remain on the singular curve in R2 defined by the equation:
5.85 − 0.00873 q 2/3 − 0.02 q 1/3 . p = q exp 3 5.85 + 0.00873 q 2/3
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
40
σu (t)
p(tf ) + 140 y(tf )
60
0.792 ≤ w ≤ 1 w = 0.7 0 ≤ w ≤ 0.628
8500
8000
20 0
7500 2000
3000
−20 0
4000
p(tf ) [mm3 ]
5
10
t [day]
(a) The Pareto front
(b) The switching function, σu (t)
15
8000 6000
0.792 ≤ w ≤ 1 w = 0.7 0 ≤ w ≤ 0.628
10
u(t)
p [mm3 ]
10
4000
5 2000
← singular control curve
0 0
5000
q [mm3 ]
10000
(c) Endothelial support, q, vs. tumor size, p
0 0
5
10
t [day] (d) Anti-angiogenic agent, u(t)
Figure 2: Tumor anti-angiogenesis: the Pareto front, and the solution curves for various weights. The singular surface allows to apply synthesis analysis which gives the following optimal bang–singular–bang control in view of the specific initial conditions in Problem (E1): for 0 ≤ t < t1 15 , u(t) = using (p(t), q(t)) , for t1 ≤ t ≤ t2 (5) 0, for t2 < t ≤ tf Since the singular control is given by a feedback expression, the computation of the bang– singular–bang control (5) is equivalent to solving the Induced Optimization Problem (cf. [32, 35]) with respect to the optimization variables (t1 , t2 , tf ). Numerically, it is more convenient to determine the optimal arc durations ξ1 = t1 , ξ2 = t2 − t1 , ξ3 = tf − t2 by the arcparametrization method; cf. [23, 32, 35]. The optimal control package NUDOCCCS developed by B¨ uskens [8, 9] allows to compute the arc durations with high precision and, moreover, provides a test of second-order sufficient conditions for the Induced Optimization Problem. Table 1 lists switching times for the specific weights displayed in Figure 2(a). Note that the switching condition (4) for the optimal control is satisfied in view of Figure 2(b), (d).
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
weights\switching times 0.792 ≤ w ≤ 1 w = 0.7 0 ≤ w ≤ 0.628
t1 1.341 1.341 1.341
t2 5.063 3.141 1.605
11
tf 9.378 8.121 7.189
Table 1: Switching times for tumor anti-angiogenesis.
4.2
Fed-batch bioreactor
We consider the problem of finding the optimal feeding rate and batch duration for a fed-batch lysine fermentation process as described by Logist et al. in [29]. The state variables and the control variable of the process are given as follows. x1 : biomass [g], x2 : substrate [g], x3 : product, lysine [g], x4 : fermenter volume [ ], u : volumetric rate of the feed stream [ /h]. Two competing objective functionals are to be maximized , in search for Pareto-optimal feeding strategy, u(t), and batch duration, tf [h] : (i) The ratio of the product formed and the process duration, i.e., the productivity : −ϕ1 (x(tf ), tf ) :=
x3 (tf ) . tf
(ii) The ratio of the product formed and the mass of the substrate added, i.e., the yield : −ϕ2 (x(tf ), tf ) :=
x3 (tf ) . 2.8 (x4 (tf ) − 5)
Note that clearly tf is free, but it is subject to the constraint 20 ≤ tf ≤ 40 . Constraints are also imposed on the fermenter volume, the feed rate and the amount of substrate to be added, respectively, as follows. 5 ≤ x4 (t) ≤ 20 , 0 ≤ u(t) ≤ 2 , 20 ≤ 2.8 (x4 (tf ) − 5) ≤ 42 . The objective functionals, constraints, and the process equations as given in [29] yield the
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
following bi-objective min subject to (E2)
optimal control problem. x3 (tf ) x3 (tf ) , − − tf 2.8 (x4 (tf ) − 5) x2 x˙ 1 = 0.125 x1 , x4 0.125 x2 x1 + 2.8 u , x˙ 2 = − 0.135 x4 x2 2 x2 + 134 0.125 x1 , x˙ 3 = −384 0.125 x4 x4 x˙ 4 = u ,
x1 (0) = 0.1 , x2 (0) = 14 , x3 (0) = 0 , x4 (0) = 5 ,
5 ≤ x4 (t) ≤ 20 ,
0 ≤ u(t) ≤ 2 , 0 ≤ tf ≤ 40 ,
12
20 ≤ 2.8 (x4 (tf ) − 5) ≤ 42 .
Scalarization of Problem (E2) can now be written as in Problem (OCPw ) as follows. min α subject to the same constraints as in (E2) and x3 (tf ) (E2w ) ∗ − β1 ≤ α , w − tf x3 (tf ) ∗ − β2 ≤ α . (1 − w) − 2.8 (x4 (tf ) − 5) Note that this is an optimal control problem with pure state constraints. The necessary optimality conditions for such problems in the form of a Minimum Principle, which can be found in [19], are more intricate than those for problems with pure control or mixed statecontrol constraints. Logist et al. [29] do not discuss the necessary optimality conditions in [19] at all. We will not go into the details of the Minimum Principle and will only sketch the main ingredients to get a better understanding of the solution structure. The Hamiltonian is defined by the direct adjoining approach in [19]: H(x, λ, µ, u) =
4
λi x˙ i + µ1 (5 − x4 ) + µ2 (x4 − 20),
(6)
i=1
where λ = (λ1 , λ2 , λ3 , λ4 ) is the adjoint variable and µ = (µ1 , µ2 ) is the multiplier associated with the state constraints 5 − x4 ≤ 0 and x4 − 20 ≤ 0. The multiplier µ satisfies the complementarity conditions µ1 (t) ≥ 0 on a boundary arc x4 (t) = 5 and µ2 (t) ≥ 0 on a boundary arc x4 (t) = 20. The adjoint equation λ˙ = −Hx is quite involved and is not written out here explicitly. Therefore, we will not give the transversality condition for λ(tf ), either. The switching function is given by σu = Hu = 2.8λ2 + λ4 . On interior arcs with 5 < x4 (t) < 20, the control minimizing the Hamiltonian is determined by the switching condition similar to (4): if σu (t) < 0 2, (7) if σu (t) > 0 u(t) = 0 , singular , if σu (t) = 0 ∀ t ∈ [t1 , t2 ] ⊂ [0, tf ]
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
13
In principle, it would be possible to compute an expression for the singular control by equating the switching function and its higher derivatives to zero. However, the calculations become so complex that we shall refrain from deriving such an expression. This is the main reason for the fact that we cannot use the arc-parametrization method for switching time optimization as in the preceding example, but rather have to employ discretization and NLP methods. On boundary arcs with either x4 (t) = 5 or x4 (t) = 20, the boundary control is determined by u = x˙ 4 = 0. In particular, the state constraints are of order one; cf., [19]. It is noteworthy that the boundary control u = 0 is on the boundary of the control set [0, 2]. In view of the minimum condition, we get the sign condition σu (t) ≥ 0
on boundary arcs.
(8)
The multiplier µ associated with the state constraints can not be computed explicitly but can be obtained from the Lagrange multiplier of the state inequality constraint of the NLP problem. We apply Algorithm 1, with N = 100 and β ∗ = (−40, −30), to generate an approximation of the Pareto front, and associated Pareto-minimum solutions for certain values of the weights as indicated in the legends, in Figure 3. In Euler discretization, we have used 2000 time nodes. This time, the Pareto front has a nonconvex part, which we can approximate well by using the Tchebychev scalarization. The “meaningful” interval of weights for this example is [0.217, 0.517]. The solution for all w ∈ [0, 0.217] is the same, corresponding to one boundary point in the Pareto front, so is the solution for all w ∈ [0.517, 1], corresponding to the other boundary point. The pure state constraint on x4 becomes active in all instances of weights; x4 attains either one or both of its lower and upper bounds during the process, i.e., x4 has one or two boundary arcs. As can be seen in Figure 3(d), the solution structure for the Pareto-optimal control is given by three types of sequences, depending on the value of the weight w: (i) 0.517 ≤ w ≤ 1: bang–bang–singular–boundary arcs (i.e., full feed–no feed–partial feed– no feed), (ii) w = 0.4: boundary–singular–boundary arcs (i.e., no feed–partial feed–no feed), (iii) w = 0.3, 0 ≤ w ≤ 0.217: boundary–singular–bang arcs (i.e., no feed–partial feed–no feed). The graphs of σu for various weights, which are depicted in Figure 3(b), match precisely the switching condition (7) on interior arcs and the sign condition (8) on boundary arcs, thus reconfirming the control structure in Figure 3(d).
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
−12
−16
0.8
σu (t)
−(yield)
1
0 ≤ w ≤ 0.217 w = 0.3 w = 0.4 0.517 ≤ w ≤ 1
−14
−18 −20
−20
−15
0.4
0 0
−10
−(productivity)
20
30
40
(b) The (scaled) switching function, σu (t)
2
15
1.5
u(t) [ /h]
20
10
0 ≤ w ≤ 0.217 w = 0.3 w = 0.4 0.517 ≤ w ≤ 1
1 0.5
5 0 0
10
t [hour]
(a) The Pareto front
x3 (t) [g/100], x4 (t) [ ]
0.6
0.2
−22 −25
14
10
20
30
40
t [hour] (c) The product, x3 (t) (scaled down by 100); the fermenter volume, x4 (t) (5 ≤ x4 (t) ≤ 20)
0 0
10
20
30
40
t [hour] (d) The feeding rate, u(t)
Figure 3: Fed-batch bioreactor: the Pareto front, and the solution curves for various weights. We apply Algorithm 1 again, with the Tchebychev scalarization replaced by the weightedsum scalarization, or equivalently the Bolza form, to Problem (E2). The result shown in Figure 4 clearly illustrates that the weighted-sum scalarization cannot recover the whole Pareto front, as expected. Moreover, the approximating points in the part of the constructed front does not look as uniformly distributed as the points in the front in Figure 3(a). Problem (E2) is also solved by Logist et al. in [29] by using the Normalized Normal Constraint scalarization technique. We improve over the results that they present on various aspects: (i) our results are obtained over a much finer discretization grid (of 2000 time nodes) and therefore they are much more accurate, (ii) we verify that the solutions we find for the control variable satisfy the necessary conditions of optimality as well as we present solutions of the key state variables.
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
15
−12
−(yield)
−14 −16 −18 −20 −22 −25
−20
−15
−10
−(productivity)
Figure 4: Fed-batch bioreactor: the Pareto front generated by the weighted-sum scalarization.
5
Conclusion
We have proposed a numerical approach for constructing an approximation of the Pareto front of nonconvex multi-objective optimal control problems. We have illustrated the working and effectiveness of our approach on two computationally challenging optimal control problems, each with two competing objectives. We have verified the necessary conditions of optimality for the discretized solutions we obtained for both problems. The fed-batch bioreactor exemplified the fact that the Bolza form, which is traditionally used (without questioning) as a way to combine conflicting objectives (or costs) with appropriate weights, may fail to represent all of the compromise, i.e. Pareto, solutions. A message of this paper therefore is that for nonconvex optimal control problems one should better use the Tchebychev scalarization instead of the Bolza form, although a differentiable re-formulation of the Tchebychev scalarization introduces a new parameter to determine and two extra terminal inequality constraints. Recall that we have defined a meaningful interval of weights to be an interval over which the Pareto solutions change qualitatively. We have worked with the meaningful interval of weights associated with the weighted Tchebychev scalarization, because it saves one significant computational effort, by having not computed the Pareto solutions outside the meaningful interval, the Pareto solutions which are the same as those found using the weights at the ends of the meaningful interval. However, the individual meaning of each weight value still remains unclear as this obscurity is an inherent property of any weight-based scalarization scheme. Therefore, it is compelling to decide on which weight, and so which solution, we are really after, by performing further optimization with an additional objective over the Pareto front, a process which is also referred to as decision making. Thus, we think that it is timely to extend the work done in [6] to nonconvex optimal control problems using the approach we have introduced in the present paper. It would be interesting to consider other scalarization techniques for multi-objective optimal control problems, especially the one recently proposed in [7] for multiobjective problems with disconnected feasible sets.
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
16
Acknowledgments The authors offer their warm thanks to the two anonymous reviewers whose comments and suggestions improved the paper.
References [1] Alt, W., Baier, R., Gerdts, M., Lempio, F.: Error bounds for Euler approximation of linear-quadratic control problems with bang–bang solutions. Numerical Algebra, Control and Optimization, 2, 547–570 (2012) [2] Alt, W., Baier, R. , Lempio, F., Gerdts, M.: Approximations of linear control problems with bang–bang solutions. Optimization, 62, 9–32 (2013) [3] Alvarez-V´azquez, L.J., Gar´ıa-Chan, N., Mart´ınez, A., V´ azquez-M´endez, M.E.: Multiobjective Pareto-optimal control: an application to wastewater management. Comput. Optim. Appl., 46, 135–157 (2010) [4] Banihashemi, N., Kaya, C.Y.: Inexact restoration for Euler discretization of boxconstrained optimal control problems, J. Optim. Theory Appl., 156, 726–760 (2013) [5] Banihashemi, N., Kaya, C.Y.: Inexact restoration and adaptive mesh refinement for constrained optimal control, J. Ind. Man. Optim., to appear. [6] Bonnel, H., Kaya, C.Y.: Optimization over the efficient set of multi-objective convex optimal control problems. J. Optim. Theory Appl., 147, 93–11 (2010). [7] Burachik, R.S., Kaya, C.Y., Rizvi, M.M.: A new scalarization technique to approximate Pareto fronts of problems with disconnected feasible sets. J. Optim. Theory Appl., DOI 10.1007/s10957-013-0346-0 (published online: 8 June 2013) [8] B¨ uskens, C.: Optimierungsmethoden und Sensitivit¨ atsanalyse f¨ ur optimale Steuerprozesse mit Steuer– und Zustands–Beschr¨ankungen. PhD Thesis, Institut f¨ ur Numerische Mathematik, Universit¨ at M¨ unster, Germany (1998) [9] B¨ uskens, C., Maurer, H.: SQP-methods for solving optimal control problems with control and state constraints: Adjoint variables, sensitivity analysis and real-time control. J. Comp. Appl. Math., 120, 85–108 (2000) [10] Craven, B.D.: Multicriteria optimal control. Asia-Pacific Journal of Operational Research, 16 53–62 (1999) [11] de Oliveira, V.A., Silva, G.N., Rojas-Medar, M.A.: A class of multiobjective control problems. Optim. Contr. Appl. Meth., 30, 77–86 (2009) [12] Dontchev, A.L., Hager, W.W.: The Euler approximation in state constrained optimal control problems. Math. Comput., 70, 173–203 (2001) [13] Dontchev, A.L., Hager, W.W., Malanowski, K.: Error bound for Euler approximation of a state and control constrained optimal control problem. Numer. Funct. Anal. Optim., 21(6), 653–682 (2000) [14] Dutta, J., Kaya, C.Y.: A new scalarization and numerical method for constructing the weak Pareto front of multi-objective optimization problems. Optimization, 60, 1091–1104 (2011)
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
17
[15] Eichfelder, G.: Adaptive Scalarization Methods in Multiobjective Optimization. Springer, Berlin, Heidelberg (2008) [16] Fourer, R., Gay, D.M., Kernighan, B.W.: AMPL: A Modeling Language for Mathematical Programming, Second Edition. Brooks/Cole Publishing Company / Cengage Learning (2003) [17] Grosset, L., Viscolani, B.: Reciprocal optimal control problems and the associated Pareto frontier. J. Optim. Theory Appl., 130, 113-123 (2006) [18] Hager, W.W., Runge-Kutta methods in optimal control and the transformed adjoint system, Numerische Mathematik, 87, 247–282 (2000) [19] Hartl, R.F., Sethi, S.P., Vickson, R.G.: A Survey of the maximum principles for optimal control problems with state constraints. SIAM Rev., 37, 181–218 (1995) [20] Jahn, J.: Vector Optimization: Theory, Applications, and Extensions, Springer-Verlag, Berlin, Heidelberg (2011) [21] Kaya, C.Y.: Inexact Restoration for Runge-Kutta discretization of optimal control problems. SIAM J. Numer. Anal., 48(4), 1492–1517 (2010) [22] Kaya, C.Y., Mart´ınez, J.M.: Euler discretization for inexact restoration and optimal control. J. Optim. Theory Appl., 134, 191–206 (2007) [23] Kaya, C.Y., Noakes, J.L.: Computational method for time-optimal switching control. J. Optim. Theory Appl., 117, 69–92 (2003) [24] Kien, B.T., Wong, N.-C., Yao, J.-C.: Necessary conditions for multiobjective optimal control problems with free end-time. SIAM J. Cont. Optim., 47(5), 2251–2274 (2010) [25] Ledzewics, U., Maurer, H., Sch¨attler, H.: Optimal and suboptimal protocols for a mathematical model for tumor anti-angiogenesis in combination with chemotherapy. Math. Biosci. Eng., 8, 307–323 (2011) [26] Ledzewics, U., Munden, J., Sch¨attler, H.: Scheduling of angiogenic inhibitors for Gomperzian and logistic tumor growth models. Discr. Contin. Dyn. Sys., Series B, 12, 415–438 (2009) [27] Ledzewics, U., Sch¨ attler, H.: Anti-angiogenic therapy in cancer treatment as an optimal control problem. SIAM J. Control Optim., 46, 1052–1079 (2007) [28] Logist, F., van Erdeghem, P.M.M., van Impe, J.F.: Efficient deterministic multiple objective optimal control of (bio)chemical processes. Chem. Eng. Sci., 64, 2527–2538 (2009) [29] Logist, F., Houska, B., Diehl, M., van Impe, J.: Fast Pareto set generation for nonlinear optimal control problems with multiple objectives. Struct. Multidisc. Optim., 42, 591-603 (2010) [30] Logist, F., Vallerio, M., Houska, B., Diehl, M., van Impe, J.: Multi-objective optimal control of chemical processes using ACADO toolkit. Comp. Chem. Eng., 37, 191–199 (2012) [31] Malanowski, K., B¨ uskens, C., Maurer, H.: Convergence of Approximations to Nonlinear Optimal Control Problems. In: Fiacco, A. V. (ed.) Mathematical Programming with Data Perturbations V, Lecture Notes in Pure and Applied Mathematics, Vol. 195, pp. 253–284. Dekker (1997)
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
18
[32] Maurer, H., B¨ uskens, C., Kim, J.-H.R., Kaya, C.Y.: Optimization methods for the verification of second-order sufficient conditions for bang–bang controls. Optim. Contr. Appl. Meth., 26, 129–156 (2005) [33] Miettinen, K.M.: Nonlinear Multiobjective Optimization, Kluwer (1999) [34] Ober-Bl¨ obaum, S., Ringkamp, M., zum Felde, G.: Solving multiobjective optimal control problems in space mission design using discrete mechanics and reference point techniques. Proceedings of the 51st IEEE Conference on Decision and Control, Dec. 10-13, Maui, Hawaii, USA, pp. 5711–5716 (2012) [35] Vossen, G.: Switching time optimization for bang-bang and singular controls. J. Optim. Theory Appl., 144, 409–429 (2006) [36] W¨achter, A., Biegler, L.T.: On the Implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Math. Progr., 106, 25–57 (2006)
Helmut Maurer†
September 3, 2013
Abstract A numerical method is proposed for constructing an approximation of the Pareto front of nonconvex multi-objective optimal control problems. First, a suitable scalarization technique is employed for the multi-objective optimal control problem. Then by using a grid of scalarization parameter values, i.e., a grid of weights, a sequence of single-objective optimal control problems are solved to obtain points which are spread over the Pareto front. The technique is illustrated on problems involving tumor anti-angiogenesis and a fed-batch bioreactor, which exhibit bang–bang, singular and boundary types of optimal control. We illustrate that the Bolza form, the traditional scalarization in optimal control, fails to represent all the compromise, i.e., Pareto optimal, solutions.
Key words: Multi-objective optimization, Optimal control, Pareto front, Numerical methods, Tumor anti-angiogenesis, Fed-batch bioreactor.
1
Introduction
In optimal control, the task is to minimize a functional, referred to as the cost, typically stated in the so-called Bolza form, tf f0 (x(t), u(t), t) dt , (1) w ψ(x(tf ), tf ) + (1 − w) 0
subject to differential and algebraic constraints, where x(t) and u(t) are the state and control variables, respectively, and w is some fixed weight selected in the interval (0, 1). The terminal time tf might be fixed or free. The terminal cost, ψ(x(tf ), tf ), is usually in conflict with the integral cost in the second term; that is why a compromise solution has to be found depending on the value of the weight w. In the multi-objective (or vector) optimization terminology, the Bolza form in (1) is the so-called weighted-sum scalarization of the multi-objective problem with the two-objective functional vector, tf f0 (x(t), u(t), t) dt , (2) ψ(x(tf ), tf ), 0
to minimize (in the Pareto sense), subject to the same constraints as those of the Bolza optimal control problem with (1). ∗
School of Information Technology and Mathematical Sciences, University of South Australia, Mawson Lakes, S.A. 5095, Australia. E-mail: [email protected] . † Institut f¨ ur Numerische und Angewandte Mathematik, Westf¨ alische Wilhelms-Universit¨ at M¨ unster, M¨ unster, Germany. E-mail: [email protected] .
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
2
Broadly speaking, the Pareto minimization of (2) is the process of finding a compromise solution, referred to as a Pareto minimum, where the value of the terminal cost or the value of the integral cost cannot be improved (i.e., reduced) further, without making the value of the other cost worse (i.e., higher). The set of all such compromise solutions form the Pareto set, or the Pareto front. While one would expect to obtain all Pareto minima by solving the single-objective Bolza problem for all w ∈ [0, 1], it turns out that this may not be possible. In other words, it does not seem to be a good idea to use the Bolza form (1) of the objective functional for nonconvex optimal control problems, as solutions found over the range of weights may not yield all of the compromise solutions. In this paper, we propose to use a scalarization different from the Bolza, or the weightedsum, form for the terminal and integral costs: we combine these costs by means of a max operator and weights, which is referred to as the weighted Tchebychev scalarization in the literature [33]. This kind of scalarization is certainly not new for general multi-objective optimization problems, see, for example, [14, 33], in the finite dimensional setting; however, to the authors’ knowledge it has not been used for optimal control problems. The weighted Tchebychev scalarization is surjective from the space of weights to the Pareto set (or Pareto front). In other words, over the range of weights, the solution of the scalarized problem yields all of the Pareto minima. Conditions of optimality and related theoretical aspects of multi-objective optimal control problems have been studied in the literature rather extensively – see for example [10, 11, 17, 24]. Moreover, effective numerical methods have been developed and illustrated for convex problems – see [3, 6] and the references therein. However, there does not seem to be many numerical studies dealing with nonconvex multi-objective optimal control problems; among existing studies we note [28, 29, 30]. In these references, first the multi-objective control problem is discretized, and then the normal boundary intersection (NBI) and normalized normal constraint (NNC) techniques are employed as a scalarization of the resulting finitedimensional problems. This process can be summarized as discretize-scalarize-then-optimize. It is worth mentioning the work in [34], where also a discretize-scalarize-then-optimize approach is adopted as in [28, 29, 30], however with the main difference that the discretization is carried out in the context of mechanics involving unconstrained dynamical equations. Optimal control problems are infinite-dimensional optimization problems, and as such the scalarization should be defined and justified for the original infinite-dimensional problem. In this paper, we choose the weighted Tchebychev scalarization for the original problem, and then consider the scalarized problem’s discretization. Therefore, our approach can be summarized as scalarize-discretize-then-optimize. In any scalarization scheme based on weights, it is hard to know the meaning of the weights used for the objective functions. It is often not even possible to know the range of values of the weights over which the solutions change qualitatively. An advantage of the weighted Tchebychev scalarization is that it is possible to determine a subinterval of the “meaningful” weights, outside of which the Pareto solutions are the same as those obtained by utilizing the weights at the boundary of the subinterval. This is a particularly useful feature (a feature which is not possessed by the weighted-sum scalarization, i.e., the Bolza form), as then the computational resources can be utilized only for finding (different) Pareto solutions corresponding to the weights in the (meaningful) subinterval mentioned above. In the NBI and NNC techniques employed in [28, 29, 30], the meaningful interval of weights would be the whole interval [0, 1]. However, the weights that are used are not associated directly with the weights of the objective functionals we consider. In the current paper, we present an algorithm for generating an approximation of the
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
3
Pareto front over the meaningful subinterval of weights. We illustrate the working of the algorithm by means of two computationally challenging problems, one involving tumor antiangiogenesis from [25] and the other a fed-batch bioreactor from [29]. For these problems, the Pareto fronts are constructed and solutions corresponding to different weights in the Pareto front are displayed. The first problem results in solutions with bang–bang and singular type of controls. The second problem involves a richer structure: the solutions, depending on the values of the weights, contain bang, singular and boundary arcs, the latter being due to the presence of (active) pure state constraints. We numerically verify that the solutions we find indeed satisfy the necessary conditions of optimality. Note that the multi-objective fed-batch bioreactor problem is also solved in [29], but over a much coarser time-grid, and the solutions found are not checked to see whether or not they satisfy the necessary conditions of optimality. We solve the multi-objective tumor anti-angiogenesis problem for the first time, in this paper. The paper is organized as follows. In Section 2, we describe the multi-objective optimal control problem, and provide definitions pertaining to scalarization. We also cite a result, which can readily be extended to the optimal control setting we have, ensuring that the Tchebychev scalarization we employ will be successful. In Section 3, we provide the algorithm. In Section 4, we implement the algorithm on the two practical optimal control problems and discuss the subsequent numerical results. In Section 5, we make concluding remarks and identify some future work.
2
Problem Statement and Preliminaries
Consider the multi-objective optimal control problem min (ϕ1 (x(tf ), tf ), . . . , ϕr (x(tf ), tf )) ˙ = f (x(t), u(t), t) , for a.e. t ∈ [0, tf ] , subject to x(t) (OCPr ) θ(x(0), x(tf ), tf ) = 0 , C(x(t), u(t), t) ≤ 0 , for a.e. t ∈ [0, tf ] , S(x(t), t) ≤ 0 , for all t ∈ [0, tf ] , where the state variable x ∈ W 1,∞ (0, tf ; IRn ), x˙ := dx/dt, and the control variable u ∈ L∞ (0, tf ; IRm ). The functions ϕi : IRn × IR+ → IR, f : IRn × IRm → IRn , θ : IRn → IRp1 , C : IRn × IRm → IRp2 , and S : IRn → IRp3 , are continuous in their arguments. In this problem, tf can be fixed or free. Here, L∞ (0, tf : IRn ) corresponds to the space of essentially bounded, measurable functions equipped with the essential supremum norm. Furthermore, W 1,∞ (0, tf ; Rn ) is the Sobolev space consisting of functions x : [0, tf ] → Rn whose first derivatives lie in L∞ . Assume that ϕi (x(tf ), tf ) ≥ 0, for all i = 1, . . . , r. Note that this assumption can easily be met by adding a large enough positive number to each objective functional. Note that Problem (OCPr ) is in general a nonsmooth problem, because it does not require differentiability of the objective functionals or the constraints. Moreover, although we have stated Problem (OCPr ) in very broad terms, it can further be generalized, for example by adding multi-point constraints, partial differential equations, etc. In other words, although Problem (OCPr ) is already in a more general form than what one usually encounters in applications, it can be further made look more general. , where tmax > 0 is some constant. Define the For technical convenience, let tf ≤ tmax f f
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
4
feasible set, X ⊂ W 1,∞ (0, tf ; IRn ) × L∞ (0, tf ; IRn ) × IR+ , such that ˙ = f (x(t), u(t), t) , for a.e. t ∈ [0, tf ]; θ(x(0), x(tf ), tf ) = 0 ; X := {(x, u, tf ) : x(t) C(x(t), u(t), t) ≤ 0 , for a.e. t ∈ [0, tf ]; S(x(t), t) ≤ 0, for all t ∈ [0, tf ]} . Define the vector of objective functionals, ϕ(x(tf ), tf ) := (ϕ1 (x(tf ), tf ), . . . , ϕr (x(tf ), tf )). The triplet (x∗ , u∗ , t∗f ) ∈ X is said to be a Pareto minimum if there exists no (x, u, tf ) ∈ X such that ϕ(x(tf ), tf ) = ϕ(x∗ (t∗f ), t∗f ) and ϕi (x(tf ), tf )) ≤ ϕi (x∗ (t∗f ), t∗f ) ,
for all i = 1, . . . , r .
On the other hand, (x∗ , u∗ , t∗f ) ∈ X is said to be a weak Pareto minimum if there exists no (x, u, tf ) ∈ X such that ϕi (x(tf ), tf )) < ϕi (x∗ (t∗f ), t∗f ) ,
for all i = 1, . . . , r .
The set of all vectors of objective functional values at the Pareto and weak Pareto minima is said to be the Pareto front (or efficient set) of Problem (OCPr ) in the r-dimensional objective value, or outcome, space. Note that the coordinates of a point in the Pareto front are simply ϕi (x∗ (t∗f ), t∗f ), i = 1, . . . , r. Obviously, when r = 2, the Pareto front is in general a curve; and when r = 3, the Pareto front is in general a surface.
2.1
Scalarization
For computing a solution of the nonconvex multi-objective optimal control problem (OCPr ), we will consider the following single-objective problem, i.e., scalarization. This scalarization is often considered for nonconvex multi-objective finite-dimensional optimization problems [33]. (Pw )
min
(x,u,tf )∈X
max{w1 ϕ1 (x(tf ), tf ), . . . , wr ϕr (x(tf ), tf )} ,
where wi , i = 1, . . . , r, are referred to as weights. The following theorem is a direct consequence of [20, Corollary 5.35]. It can also be concluded by generalizing the proofs given in a finite-dimensional setting in [33, Theorems 3.4.2 and 3.4.5] to the infinite dimensional setting in a straightforward manner. Theorem 1 The triplet (x∗ , u∗ , t∗f ) is a weak Pareto minimum of (OCPr ) if, and only if, (x∗ , u∗ , t∗f ) is a solution of (Pw ) for some w1 , . . . , wr > 0. An ideal cost ϕ∗i , i = 1, . . . , p, associated with Problem (Pw ) is the optimal value of the optimal control problem, min ϕi (x(tf ), tf ) . (Pi ) (x,u,tf )∈X
Let (x, u, tf ) be a minimizer of the single-objective problem (Pi ) above. Then we also define ϕj := ϕj (x(tf ), tf ), for j = i and j = 1, . . . , r. Problem (Pw ) is referred to as the weighted Tchebychev problem (or Tchebychev scalarization) because of the weighted Tchebychev norm, maxi |wi ϕi (x(tf ), tf )| = maxi wi ϕi (x(tf ), tf ), appearing in the objective. In the case when ϕ∗i is negative, one can simply add a large enough positive number to the ith objective, to make the objective positive. In general, it is useful to add a positive number to each objective in order to obtain an even spread of the Pareto points approximating the
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
5
Pareto front – see, for example, [14], for further discussion and geometric illustration. For this purpose, we define next the so-called utopian objective values. A utopian objective vector β ∗ associated with Problem (OCPr ) consists of the components given as βi∗ = ϕ∗i − εi where εi > 0 for all i = 1, . . . , p. Problem (Pw ) can be equivalently written as βi∗
min
(x,u,tf )∈X
max{w1 (ϕ1 (x(tf ), tf ) − β1∗ ), . . . , wr (ϕr (x(tf ), tf ) − βr∗ )} .
In the case when the objective functionals and the constraints in Problem (OCPr ) are differentiable in their arguments, Problem (Pw ), or its equivalent above, is still a non-smooth optimization problem, because of the max operator in the objective. However, it is well-known that Problem (Pw ) can be re-formulated in the following (smooth) form. min α α≥0 (x,u,t )∈X f subject to w1 (ϕ1 (x(tf ), tf ) − β1∗ ) ≤ α , (OCPw ) .. . wr (ϕr (x(tf ), tf ) − βr∗ ) ≤ α . Problem (OCPw ) is referred to as goal attainment method [33], as well as Pascoletti-Serafini scalarization [15]. We will solve Problem (OCPw ) in an algorithm we present in the next section, for the two examples we want to study. Note that the “popular” weighted-sum scalarization, given below, fails to generate the nonconvex parts of a Pareto front. (Pws )
min
(x,u,tf )∈X
r
wi ϕi (x(tf ), tf ) .
i=1
We will illustrate this in the fed-batch bioreactor problem in Section 4.2. It should be noted that, for r = 2, the above weighted-sum scalarization can be shown to be equivalent to the Bolza form in (1), after defining a new state variable, as in the case of the tumor antiangiogenesis problem in Section 4.1.
3
An Algorithm to Generate the Pareto Front
In this section, we provide an algorithm which generates an approximation of the Pareto front of Problem (OCPr ) with two objectives. Note that by choosing w1 = w, and w2 = 1 − w, where w ∈ [0, 1], one can simply consider the single weight w in the scalarized problem (OCPw ). Solution of Problem (OCPw ) using equal increments of the value of w in [0, 1] can generate a reasonable approximation of the Pareto front, if one chooses a utopia point (β1∗ , β2∗ ) far enough from the front (see [14]). With the Tchebychev scalarization, the weight w would usually take values over a (smaller) subinterval of [0, 1]. In other words, a restricted range of values of w, i.e., 0 < w0 < w < wf < 1 would usually suffice for the generation of the front. Figure 3 illustrates the geometry with which one can compute the subinterval end-points, w0 and wf . Note that the points (ϕ∗1 , ϕ2 ) and (ϕ1 , ϕ∗2 ) represent the boundary of the Pareto front, where the equations of the “rays” emanating from the utopia point (β1∗ , β2∗ ) and passing through the boundary points are given as in the figure. Also note that these equations are obtained by virtue of the weighted
6
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
ϕ2
wf (ϕ1 − β1∗ ) = (1 − wf ) (ϕ2 − β2∗ ) (ϕ∗1 , ϕ2 ) Pareto front w0 (ϕ1 − β1∗ ) = (1 − w0 ) (ϕ2 − β2∗ ) (ϕ1 , ϕ∗2 ) ϕ1
(β1∗ , β2∗ )
Figure 1: Determination of the boundary weights w0 and wf .
Tchebychev (or max) norm. By substituting the boundary values of the Pareto curve into the respective equations, and solving each equation for w0 and wf one simply gets w0 =
(ϕ∗2 − β2∗ ) (ϕ1 − β1∗ ) + (ϕ∗2 − β2∗ )
and
wf =
(ϕ2 − β2∗ ) . (ϕ∗1 − β1∗ ) + (ϕ2 − β2∗ )
(3)
Although the algorithm below is described for the bi-objective case (i.e., when r = 2), it can be generalized to problems with more than two objectives, by re-defining the boundary of the space of weights and generating a grid of weights on a segment of a hyperplane rather than a subinterval. However, care needs be taken in doing this: for example for r = 3, it is well known that the boundary of the Pareto front (which now typically is a closed curve rather than just two points) may very well lie outside the triangle formed by the three points, (ϕ∗1 , ϕ2 , ϕ3 ), (ϕ1 , ϕ∗2 , ϕ3 ) and (ϕ1 , ϕ2 , ϕ∗3 ). While this kind of generalization is quite important, it is outside the scope of the present paper. Algorithm 1 Step 0.0 (Initialization) Choose the utopia parameters, ε1 , ε2 > 0. Set the number of discrete (approximating) points, (N + 1), in the Pareto front. Set k := 1. Step 0.1 (Boundary of the front) 1
1
1
2
2
2
(a) Find (x1 , u1 , tf ) that solves Problem (P1 ). Note that ϕ∗1 = ϕ1 (x1 (tf ), tf ) and 1 1 ϕ2 = ϕ2 (x1 (tf ), tf ). Mark a boundary point in the Pareto front: ϕN := (ϕ∗1 , ϕ2 ) . (b) Find (x2 , u2 , tf ) that solves Problem (P2 ). Note that ϕ1 = ϕ1 (x2 (tf ), tf ) and 2 2 ϕ∗2 = ϕ2 (x2 (tf ), tf ). Mark a boundary point in the Pareto front: ϕ0 := (ϕ1 , ϕ∗2 ) .
Step 0.2 (Utopia point) Set β ∗ := (β1∗ , β2∗ ) with βi∗ := ϕ∗i − εi , i = 1, 2. Step 0.3 (Range of weights) Determine the range of (meaningful) weights, i.e., the subinterval [w0 , wf ], using (3). Set the increment, ∆w := (wf − w0 )/N .
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
7
Step k.1 (Current weights) Set w := w0 + k ∆w. Set w1 := w and w2 := 1 − w. Step k.2 (A Pareto minimum) Find (x∗ , u∗ , t∗f ) that solves Problem (OCPw ). Assign a point in the Pareto front: ϕk := (ϕ1 (x∗ (t∗f ), t∗f ), ϕ2 (x∗ (t∗f ), t∗f )). Step k.3 (Stopping criterion) If k = N then STOP. Otherwise, set k := k + 1 , and go to Step k.1. In Step k.1 of Algorithm 1, the weights are chosen equidistantly. One should note that for a uniform distribution of the points in the value space approximating the Pareto front, these weights can as well be chosen in an adaptive fashion as discussed by Eichfelder in [15].
4
Numerical Examples
In this section, we illustrate the working of the algorithm on two numerical example problems, one involving tumor anti-angiogenesis in Section 4.1 and the other a fed-batch bioreactor in Section 4.2. It is common practice to approximate the state and control variables over a partition of the time horizon in optimal control problems, a process referred to as discretization, so that one can solve a (large-scale) finite-dimensional optimization problem, instead of an infinitedimensional one, to obtain an approximate solution – see, for example, [4, 5, 12, 13, 18, 21, 22, 31]. This overall process is referred to as a direct method, or a discretize-then-optimize approach. Convergence theory for the discretize-then-optimize approach and its applications have been studied in the literature extensively; see the above-listed references, which in general concern Runge-Kutta (and, as a particular case, Euler) discretization of state- and controlconstrained optimal control problems. In [12], where Euler discretization is considered, for example, under the rather general assumptions that the control variable is continuous and second-order sufficient conditions hold, the discrete (i.e., approximate) solutions of the state, control and adjoint variables converge to the continuous-time solutions (the state and costate variables with respect to the W 1,∞ -norm and the control variable with respect to the L∞ norm). Higher-order Runge-Kutta discretization schemes, which require higher-order differentiability of the solution, in general provide higher-order convergence rates [18]. In the fed-batch bioreactor example, we discretize the problem over a sufficiently fine grid. The solution for the control variable is discontinuous, incorporating bang–bang, singular and boundary arcs. So far, there exists no formal proof of convergence that the discretized solution for this example converges to the infinite-dimensional solution in topology of the L1 space. Such a convergence result has been obtained by Alt et al. [1, 2] only in the case of a linear dynamics and a bang–bang solution. A convergence result for more general problems involving discontinuous solutions has not been established yet in the literature. Nevertheless, our discretization approach clearly shows the control structure of the fed-batch bioreactor example. In the tumor anti-angiogenesis example, the solution also involves bang–bang and singular controls. However, since we can derive a feedback expression for the candidate singular optimal control, we parameterize the problem with respect to the (unknown) times when the bang–bang and singular controls switch, and so we reduce the scalarized problem into what is called an Induced Problem, which is explained in some detail in Section 4.1. We still employ a discretization scheme for computing the bang-bang and singular arcs precisely enough.
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
8
In Steps 0.1 and k.2 of the algorithm for the example problems that we study, a discretized version of Problem (OCPw ) is solved by using Ipopt, version 3.2.4s, a popular optimization software based on an interior point method; see [36]. We use AMPL [16] as an optimization modelling language, which employs Ipopt as a solver. Note that, by Theorem 1, there is a bijection between the space of weights and the space of Pareto minima. So, as the grid of weights is refined, one would expect to find more points in the Pareto front and obtain a better approximation, ultimately converging to the true Pareto front. Recall that here we take the scalarize-discretize-then-optimize approach, meaning that we first scalarize the multi-objective problem in the infinite-dimensional space, and then discretize before we apply a finite-dimensional optimization method to find a discrete approximate solution of the scalarized problem. By the existing theory of discretization we outlined above, under certain assumptions, one can guarantee that (or at least check the necessary and/or sufficient optimality conditions a posteriori to observe that) the discrete approximate solution converges to a solution of the continuous-time scalarization of the original problem, yielding a Pareto minimum of the original problem. On the other hand, if one takes the existing discretize-scalarize-then-optimize approach, then scalarization is carried out for the discretization of the original multi-objective problem. To the authors’ knowledge, no theory is available for guaranteeing convergence of a solution of the discretized multi-objective control problem to a solution of the continuous-time multiobjective control problem.
4.1
Optimal protocols for tumor anti-angiogenesis
Let us cite from the Introduction in the article of Ledzewicz, Maurer and Sch¨attler [25]: “Tumor anti-angiogenesis is an indirect cancer treatment approach with the aim to limit the tumor’s growth, possibly even shrink the tumor, by depriving it of the vasculature it needs for a steady supply with nutrients.” In a seminal paper, Ledzewicz and Sch¨attler [27] formulated the anti-angiogenic cancer treatment as an optimal control problem and identified the optimal control structure as a bang–singular–bang control. Here, we consider the dynamical model for anti-angiogenic monotherapy in [25] which is slightly different from the one in [27]. The state and control variables of the problem are as follows. p : primary tumor volume [mm3 ], q : carrying capacity, or endothelial support [mm3 ], u : anti-angiogenic agent. We introduce two competing objective functionals to minimize : (i) The final tumor volume : p(tf ) , (ii) The final tumor volume, plus a factor of the total amount of anti-angiogenic agent administered : tf u(t) dt . p(tf ) + 140 0
The factor 140 in the second functional has been employed to guarantee the same optimal control structure for all weights of interest. Note that the terminal time tf is free. It may t at first seem that it is “more natural” to consider 0 f u(t) dt alone as the second objective functional; however, the minimization of this integral cost (for finding a boundary point of
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
9
the front) yields the trivial solution u(t) = 0 and tf = 0 (tf is free), which is not realistic. Therefore we consider a combination of p(tf ) and the integral cost with a large positive coefficient. Incorporating the dynamical equations and constraints as described following bi-objective optimal control problem: min (p(tf ), p(tf ) + 140 y(tf )) p subject to p˙ = −0.084 p ln q , (E1) q˙ = 5.85 q 2/3 − 0.00873 q 4/3 − 0.02 q − 0.15 q u , y˙ = u , y(tf ) ≤ 45 , 0 ≤ u(t) ≤ 15 .
in [25], we get the
p(0) = 8000 , q(0) = 10000 , y(0) = 0 ,
The scalarization of Problem (E1) in the form of problem (OCPw ) in Section 2 yields the following optimal control problem. min α subject to the same constraints as in (E1) and (E1w ) w (p(tf ) − β1∗ ) ≤ α , (1 − w) (p(tf ) + 140 y(tf ) − β2∗ ) ≤ α . Note that, in the above problem, ϕ2 ≥ ϕ1 . This has a geometric implication on the Pareto front and the weights that are associated in constructing it: the front simply lies above the line ϕ2 = ϕ1 in the ϕ1 -ϕ2 plane, and the meaningful interval of weights is, for β ∗ = (0, 0), a subinterval of [0.5, 1]. Algorithm 1, with N = 50 and β ∗ = (0, 0), results in the approximation of the Pareto front shown in Figure 2(a). The “meaningful” interval of weights is [0.628, 0.792], which is a subset of [0.5, 1] as expected. The solution for all w ∈ [0, 0.628] is the same, corresponding to one boundary point in the Pareto front, so is the solution for all w ∈ [0.792, 1], corresponding to the other boundary point. The switching function σu (t) = −0.15 λq (t) q(t) + λy (t) , where λq and λy are the adjoint (or costate) variables associated with the differential equations for the state variables q and y, determines the minimizing control according to if σu (t) < 0 15 , (4) if σu (t) > 0 u(t) = 0 , singular , if σu (t) = 0 ∀ t ∈ [t1 , t2 ] ⊂ [0, tf ] ¨u (t) = 0 ∀ t ∈ [t1 , t2 ], A singular control can be computed from the relations σu (t) = σ˙ u (t) = σ since any singular arc has order one. In [25, 26], the following feedback expression was obtained for a singular control:
5.85 − 0.00873 q 2/3 1 5.85 + 0.00873 q 2/3 + 3 (0.084) − 0.02 . using (p, q) = 0.02 q 1/3 5.85 − 0.00873 q 2/3 Moreover, all singular trajectories remain on the singular curve in R2 defined by the equation:
5.85 − 0.00873 q 2/3 − 0.02 q 1/3 . p = q exp 3 5.85 + 0.00873 q 2/3
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
40
σu (t)
p(tf ) + 140 y(tf )
60
0.792 ≤ w ≤ 1 w = 0.7 0 ≤ w ≤ 0.628
8500
8000
20 0
7500 2000
3000
−20 0
4000
p(tf ) [mm3 ]
5
10
t [day]
(a) The Pareto front
(b) The switching function, σu (t)
15
8000 6000
0.792 ≤ w ≤ 1 w = 0.7 0 ≤ w ≤ 0.628
10
u(t)
p [mm3 ]
10
4000
5 2000
← singular control curve
0 0
5000
q [mm3 ]
10000
(c) Endothelial support, q, vs. tumor size, p
0 0
5
10
t [day] (d) Anti-angiogenic agent, u(t)
Figure 2: Tumor anti-angiogenesis: the Pareto front, and the solution curves for various weights. The singular surface allows to apply synthesis analysis which gives the following optimal bang–singular–bang control in view of the specific initial conditions in Problem (E1): for 0 ≤ t < t1 15 , u(t) = using (p(t), q(t)) , for t1 ≤ t ≤ t2 (5) 0, for t2 < t ≤ tf Since the singular control is given by a feedback expression, the computation of the bang– singular–bang control (5) is equivalent to solving the Induced Optimization Problem (cf. [32, 35]) with respect to the optimization variables (t1 , t2 , tf ). Numerically, it is more convenient to determine the optimal arc durations ξ1 = t1 , ξ2 = t2 − t1 , ξ3 = tf − t2 by the arcparametrization method; cf. [23, 32, 35]. The optimal control package NUDOCCCS developed by B¨ uskens [8, 9] allows to compute the arc durations with high precision and, moreover, provides a test of second-order sufficient conditions for the Induced Optimization Problem. Table 1 lists switching times for the specific weights displayed in Figure 2(a). Note that the switching condition (4) for the optimal control is satisfied in view of Figure 2(b), (d).
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
weights\switching times 0.792 ≤ w ≤ 1 w = 0.7 0 ≤ w ≤ 0.628
t1 1.341 1.341 1.341
t2 5.063 3.141 1.605
11
tf 9.378 8.121 7.189
Table 1: Switching times for tumor anti-angiogenesis.
4.2
Fed-batch bioreactor
We consider the problem of finding the optimal feeding rate and batch duration for a fed-batch lysine fermentation process as described by Logist et al. in [29]. The state variables and the control variable of the process are given as follows. x1 : biomass [g], x2 : substrate [g], x3 : product, lysine [g], x4 : fermenter volume [ ], u : volumetric rate of the feed stream [ /h]. Two competing objective functionals are to be maximized , in search for Pareto-optimal feeding strategy, u(t), and batch duration, tf [h] : (i) The ratio of the product formed and the process duration, i.e., the productivity : −ϕ1 (x(tf ), tf ) :=
x3 (tf ) . tf
(ii) The ratio of the product formed and the mass of the substrate added, i.e., the yield : −ϕ2 (x(tf ), tf ) :=
x3 (tf ) . 2.8 (x4 (tf ) − 5)
Note that clearly tf is free, but it is subject to the constraint 20 ≤ tf ≤ 40 . Constraints are also imposed on the fermenter volume, the feed rate and the amount of substrate to be added, respectively, as follows. 5 ≤ x4 (t) ≤ 20 , 0 ≤ u(t) ≤ 2 , 20 ≤ 2.8 (x4 (tf ) − 5) ≤ 42 . The objective functionals, constraints, and the process equations as given in [29] yield the
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
following bi-objective min subject to (E2)
optimal control problem. x3 (tf ) x3 (tf ) , − − tf 2.8 (x4 (tf ) − 5) x2 x˙ 1 = 0.125 x1 , x4 0.125 x2 x1 + 2.8 u , x˙ 2 = − 0.135 x4 x2 2 x2 + 134 0.125 x1 , x˙ 3 = −384 0.125 x4 x4 x˙ 4 = u ,
x1 (0) = 0.1 , x2 (0) = 14 , x3 (0) = 0 , x4 (0) = 5 ,
5 ≤ x4 (t) ≤ 20 ,
0 ≤ u(t) ≤ 2 , 0 ≤ tf ≤ 40 ,
12
20 ≤ 2.8 (x4 (tf ) − 5) ≤ 42 .
Scalarization of Problem (E2) can now be written as in Problem (OCPw ) as follows. min α subject to the same constraints as in (E2) and x3 (tf ) (E2w ) ∗ − β1 ≤ α , w − tf x3 (tf ) ∗ − β2 ≤ α . (1 − w) − 2.8 (x4 (tf ) − 5) Note that this is an optimal control problem with pure state constraints. The necessary optimality conditions for such problems in the form of a Minimum Principle, which can be found in [19], are more intricate than those for problems with pure control or mixed statecontrol constraints. Logist et al. [29] do not discuss the necessary optimality conditions in [19] at all. We will not go into the details of the Minimum Principle and will only sketch the main ingredients to get a better understanding of the solution structure. The Hamiltonian is defined by the direct adjoining approach in [19]: H(x, λ, µ, u) =
4
λi x˙ i + µ1 (5 − x4 ) + µ2 (x4 − 20),
(6)
i=1
where λ = (λ1 , λ2 , λ3 , λ4 ) is the adjoint variable and µ = (µ1 , µ2 ) is the multiplier associated with the state constraints 5 − x4 ≤ 0 and x4 − 20 ≤ 0. The multiplier µ satisfies the complementarity conditions µ1 (t) ≥ 0 on a boundary arc x4 (t) = 5 and µ2 (t) ≥ 0 on a boundary arc x4 (t) = 20. The adjoint equation λ˙ = −Hx is quite involved and is not written out here explicitly. Therefore, we will not give the transversality condition for λ(tf ), either. The switching function is given by σu = Hu = 2.8λ2 + λ4 . On interior arcs with 5 < x4 (t) < 20, the control minimizing the Hamiltonian is determined by the switching condition similar to (4): if σu (t) < 0 2, (7) if σu (t) > 0 u(t) = 0 , singular , if σu (t) = 0 ∀ t ∈ [t1 , t2 ] ⊂ [0, tf ]
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
13
In principle, it would be possible to compute an expression for the singular control by equating the switching function and its higher derivatives to zero. However, the calculations become so complex that we shall refrain from deriving such an expression. This is the main reason for the fact that we cannot use the arc-parametrization method for switching time optimization as in the preceding example, but rather have to employ discretization and NLP methods. On boundary arcs with either x4 (t) = 5 or x4 (t) = 20, the boundary control is determined by u = x˙ 4 = 0. In particular, the state constraints are of order one; cf., [19]. It is noteworthy that the boundary control u = 0 is on the boundary of the control set [0, 2]. In view of the minimum condition, we get the sign condition σu (t) ≥ 0
on boundary arcs.
(8)
The multiplier µ associated with the state constraints can not be computed explicitly but can be obtained from the Lagrange multiplier of the state inequality constraint of the NLP problem. We apply Algorithm 1, with N = 100 and β ∗ = (−40, −30), to generate an approximation of the Pareto front, and associated Pareto-minimum solutions for certain values of the weights as indicated in the legends, in Figure 3. In Euler discretization, we have used 2000 time nodes. This time, the Pareto front has a nonconvex part, which we can approximate well by using the Tchebychev scalarization. The “meaningful” interval of weights for this example is [0.217, 0.517]. The solution for all w ∈ [0, 0.217] is the same, corresponding to one boundary point in the Pareto front, so is the solution for all w ∈ [0.517, 1], corresponding to the other boundary point. The pure state constraint on x4 becomes active in all instances of weights; x4 attains either one or both of its lower and upper bounds during the process, i.e., x4 has one or two boundary arcs. As can be seen in Figure 3(d), the solution structure for the Pareto-optimal control is given by three types of sequences, depending on the value of the weight w: (i) 0.517 ≤ w ≤ 1: bang–bang–singular–boundary arcs (i.e., full feed–no feed–partial feed– no feed), (ii) w = 0.4: boundary–singular–boundary arcs (i.e., no feed–partial feed–no feed), (iii) w = 0.3, 0 ≤ w ≤ 0.217: boundary–singular–bang arcs (i.e., no feed–partial feed–no feed). The graphs of σu for various weights, which are depicted in Figure 3(b), match precisely the switching condition (7) on interior arcs and the sign condition (8) on boundary arcs, thus reconfirming the control structure in Figure 3(d).
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
−12
−16
0.8
σu (t)
−(yield)
1
0 ≤ w ≤ 0.217 w = 0.3 w = 0.4 0.517 ≤ w ≤ 1
−14
−18 −20
−20
−15
0.4
0 0
−10
−(productivity)
20
30
40
(b) The (scaled) switching function, σu (t)
2
15
1.5
u(t) [ /h]
20
10
0 ≤ w ≤ 0.217 w = 0.3 w = 0.4 0.517 ≤ w ≤ 1
1 0.5
5 0 0
10
t [hour]
(a) The Pareto front
x3 (t) [g/100], x4 (t) [ ]
0.6
0.2
−22 −25
14
10
20
30
40
t [hour] (c) The product, x3 (t) (scaled down by 100); the fermenter volume, x4 (t) (5 ≤ x4 (t) ≤ 20)
0 0
10
20
30
40
t [hour] (d) The feeding rate, u(t)
Figure 3: Fed-batch bioreactor: the Pareto front, and the solution curves for various weights. We apply Algorithm 1 again, with the Tchebychev scalarization replaced by the weightedsum scalarization, or equivalently the Bolza form, to Problem (E2). The result shown in Figure 4 clearly illustrates that the weighted-sum scalarization cannot recover the whole Pareto front, as expected. Moreover, the approximating points in the part of the constructed front does not look as uniformly distributed as the points in the front in Figure 3(a). Problem (E2) is also solved by Logist et al. in [29] by using the Normalized Normal Constraint scalarization technique. We improve over the results that they present on various aspects: (i) our results are obtained over a much finer discretization grid (of 2000 time nodes) and therefore they are much more accurate, (ii) we verify that the solutions we find for the control variable satisfy the necessary conditions of optimality as well as we present solutions of the key state variables.
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
15
−12
−(yield)
−14 −16 −18 −20 −22 −25
−20
−15
−10
−(productivity)
Figure 4: Fed-batch bioreactor: the Pareto front generated by the weighted-sum scalarization.
5
Conclusion
We have proposed a numerical approach for constructing an approximation of the Pareto front of nonconvex multi-objective optimal control problems. We have illustrated the working and effectiveness of our approach on two computationally challenging optimal control problems, each with two competing objectives. We have verified the necessary conditions of optimality for the discretized solutions we obtained for both problems. The fed-batch bioreactor exemplified the fact that the Bolza form, which is traditionally used (without questioning) as a way to combine conflicting objectives (or costs) with appropriate weights, may fail to represent all of the compromise, i.e. Pareto, solutions. A message of this paper therefore is that for nonconvex optimal control problems one should better use the Tchebychev scalarization instead of the Bolza form, although a differentiable re-formulation of the Tchebychev scalarization introduces a new parameter to determine and two extra terminal inequality constraints. Recall that we have defined a meaningful interval of weights to be an interval over which the Pareto solutions change qualitatively. We have worked with the meaningful interval of weights associated with the weighted Tchebychev scalarization, because it saves one significant computational effort, by having not computed the Pareto solutions outside the meaningful interval, the Pareto solutions which are the same as those found using the weights at the ends of the meaningful interval. However, the individual meaning of each weight value still remains unclear as this obscurity is an inherent property of any weight-based scalarization scheme. Therefore, it is compelling to decide on which weight, and so which solution, we are really after, by performing further optimization with an additional objective over the Pareto front, a process which is also referred to as decision making. Thus, we think that it is timely to extend the work done in [6] to nonconvex optimal control problems using the approach we have introduced in the present paper. It would be interesting to consider other scalarization techniques for multi-objective optimal control problems, especially the one recently proposed in [7] for multiobjective problems with disconnected feasible sets.
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
16
Acknowledgments The authors offer their warm thanks to the two anonymous reviewers whose comments and suggestions improved the paper.
References [1] Alt, W., Baier, R., Gerdts, M., Lempio, F.: Error bounds for Euler approximation of linear-quadratic control problems with bang–bang solutions. Numerical Algebra, Control and Optimization, 2, 547–570 (2012) [2] Alt, W., Baier, R. , Lempio, F., Gerdts, M.: Approximations of linear control problems with bang–bang solutions. Optimization, 62, 9–32 (2013) [3] Alvarez-V´azquez, L.J., Gar´ıa-Chan, N., Mart´ınez, A., V´ azquez-M´endez, M.E.: Multiobjective Pareto-optimal control: an application to wastewater management. Comput. Optim. Appl., 46, 135–157 (2010) [4] Banihashemi, N., Kaya, C.Y.: Inexact restoration for Euler discretization of boxconstrained optimal control problems, J. Optim. Theory Appl., 156, 726–760 (2013) [5] Banihashemi, N., Kaya, C.Y.: Inexact restoration and adaptive mesh refinement for constrained optimal control, J. Ind. Man. Optim., to appear. [6] Bonnel, H., Kaya, C.Y.: Optimization over the efficient set of multi-objective convex optimal control problems. J. Optim. Theory Appl., 147, 93–11 (2010). [7] Burachik, R.S., Kaya, C.Y., Rizvi, M.M.: A new scalarization technique to approximate Pareto fronts of problems with disconnected feasible sets. J. Optim. Theory Appl., DOI 10.1007/s10957-013-0346-0 (published online: 8 June 2013) [8] B¨ uskens, C.: Optimierungsmethoden und Sensitivit¨ atsanalyse f¨ ur optimale Steuerprozesse mit Steuer– und Zustands–Beschr¨ankungen. PhD Thesis, Institut f¨ ur Numerische Mathematik, Universit¨ at M¨ unster, Germany (1998) [9] B¨ uskens, C., Maurer, H.: SQP-methods for solving optimal control problems with control and state constraints: Adjoint variables, sensitivity analysis and real-time control. J. Comp. Appl. Math., 120, 85–108 (2000) [10] Craven, B.D.: Multicriteria optimal control. Asia-Pacific Journal of Operational Research, 16 53–62 (1999) [11] de Oliveira, V.A., Silva, G.N., Rojas-Medar, M.A.: A class of multiobjective control problems. Optim. Contr. Appl. Meth., 30, 77–86 (2009) [12] Dontchev, A.L., Hager, W.W.: The Euler approximation in state constrained optimal control problems. Math. Comput., 70, 173–203 (2001) [13] Dontchev, A.L., Hager, W.W., Malanowski, K.: Error bound for Euler approximation of a state and control constrained optimal control problem. Numer. Funct. Anal. Optim., 21(6), 653–682 (2000) [14] Dutta, J., Kaya, C.Y.: A new scalarization and numerical method for constructing the weak Pareto front of multi-objective optimization problems. Optimization, 60, 1091–1104 (2011)
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
17
[15] Eichfelder, G.: Adaptive Scalarization Methods in Multiobjective Optimization. Springer, Berlin, Heidelberg (2008) [16] Fourer, R., Gay, D.M., Kernighan, B.W.: AMPL: A Modeling Language for Mathematical Programming, Second Edition. Brooks/Cole Publishing Company / Cengage Learning (2003) [17] Grosset, L., Viscolani, B.: Reciprocal optimal control problems and the associated Pareto frontier. J. Optim. Theory Appl., 130, 113-123 (2006) [18] Hager, W.W., Runge-Kutta methods in optimal control and the transformed adjoint system, Numerische Mathematik, 87, 247–282 (2000) [19] Hartl, R.F., Sethi, S.P., Vickson, R.G.: A Survey of the maximum principles for optimal control problems with state constraints. SIAM Rev., 37, 181–218 (1995) [20] Jahn, J.: Vector Optimization: Theory, Applications, and Extensions, Springer-Verlag, Berlin, Heidelberg (2011) [21] Kaya, C.Y.: Inexact Restoration for Runge-Kutta discretization of optimal control problems. SIAM J. Numer. Anal., 48(4), 1492–1517 (2010) [22] Kaya, C.Y., Mart´ınez, J.M.: Euler discretization for inexact restoration and optimal control. J. Optim. Theory Appl., 134, 191–206 (2007) [23] Kaya, C.Y., Noakes, J.L.: Computational method for time-optimal switching control. J. Optim. Theory Appl., 117, 69–92 (2003) [24] Kien, B.T., Wong, N.-C., Yao, J.-C.: Necessary conditions for multiobjective optimal control problems with free end-time. SIAM J. Cont. Optim., 47(5), 2251–2274 (2010) [25] Ledzewics, U., Maurer, H., Sch¨attler, H.: Optimal and suboptimal protocols for a mathematical model for tumor anti-angiogenesis in combination with chemotherapy. Math. Biosci. Eng., 8, 307–323 (2011) [26] Ledzewics, U., Munden, J., Sch¨attler, H.: Scheduling of angiogenic inhibitors for Gomperzian and logistic tumor growth models. Discr. Contin. Dyn. Sys., Series B, 12, 415–438 (2009) [27] Ledzewics, U., Sch¨ attler, H.: Anti-angiogenic therapy in cancer treatment as an optimal control problem. SIAM J. Control Optim., 46, 1052–1079 (2007) [28] Logist, F., van Erdeghem, P.M.M., van Impe, J.F.: Efficient deterministic multiple objective optimal control of (bio)chemical processes. Chem. Eng. Sci., 64, 2527–2538 (2009) [29] Logist, F., Houska, B., Diehl, M., van Impe, J.: Fast Pareto set generation for nonlinear optimal control problems with multiple objectives. Struct. Multidisc. Optim., 42, 591-603 (2010) [30] Logist, F., Vallerio, M., Houska, B., Diehl, M., van Impe, J.: Multi-objective optimal control of chemical processes using ACADO toolkit. Comp. Chem. Eng., 37, 191–199 (2012) [31] Malanowski, K., B¨ uskens, C., Maurer, H.: Convergence of Approximations to Nonlinear Optimal Control Problems. In: Fiacco, A. V. (ed.) Mathematical Programming with Data Perturbations V, Lecture Notes in Pure and Applied Mathematics, Vol. 195, pp. 253–284. Dekker (1997)
A Numerical Method for Nonconvex Multi-objective Optimal Control Problems by C. Y. Kaya and H. Maurer
18
[32] Maurer, H., B¨ uskens, C., Kim, J.-H.R., Kaya, C.Y.: Optimization methods for the verification of second-order sufficient conditions for bang–bang controls. Optim. Contr. Appl. Meth., 26, 129–156 (2005) [33] Miettinen, K.M.: Nonlinear Multiobjective Optimization, Kluwer (1999) [34] Ober-Bl¨ obaum, S., Ringkamp, M., zum Felde, G.: Solving multiobjective optimal control problems in space mission design using discrete mechanics and reference point techniques. Proceedings of the 51st IEEE Conference on Decision and Control, Dec. 10-13, Maui, Hawaii, USA, pp. 5711–5716 (2012) [35] Vossen, G.: Switching time optimization for bang-bang and singular controls. J. Optim. Theory Appl., 144, 409–429 (2006) [36] W¨achter, A., Biegler, L.T.: On the Implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Math. Progr., 106, 25–57 (2006)
