A Global Approach to the Optimal Control of System
Dynamics Models
A. Fiigenschuh!, I. Vierhaus?
1 Helmut Schmidt University / University of the Federal Armed Forces Hamburg, Germany
2 Zuse Institute Berlin, Germany
Abstract:
‘The System Dynamics (SD) methodology is a framework for modeling and simulating the dynamic behavior
the occurrence of feedback
be the system are usually
simple systems
can show a nonintuitive, unpredictable behavior over time. Controlling a dynamical system means to
specify potential interventions from outside that should keep the system on the desired track, and to
define an evaluation
of socioeconomic systems. Characteristic for the description of such systems
loops together with stocks and flows. The mathematical equations that de
ordinary differential equations and nonlinear algebraic constraints. ‘Therefore seemingly
hema to compare different controls among each other, so that a “best” control can
be defined in a meaningful way. The central question is how to compute such globally optimal control
for a given SD model, that allows the transition of the system into a desired state with minimum effort,
We propose a nonlinear (MINLP) of the System Dynamics
Optimization (SDO) sae, MINED protleus cat be solved by linear progtamamning based beahek aad:
bound approach. We demonstrate that standard MINLP solvers are not able to solve SDO problem. To
overcome this obstacle, we introduce a special-tailored bound propagation method. We apply our new
method to a predator-prey model with additional hunting activity as control, and to a mini-world model
with the consumption level as control. Numerical results for these test cases are presented
Keywords:
Global Optimal Control; Mixed-Integer Nonlinear Optimization; Bounds
System Dynami
ing.
1 Introduction: From SD to MINLP
SD models des
loops. The relation between those
functional relations, logical relations (such as i
is individually well understood, the interplay of several of these relations can show a surp:
ibe the behavior of a
stem that consists of several interrelated stocks, flows and feedback
usually determined by ordinary differential equations, nonlinear
-then-clse), or tabular data. Even if each of these relations
ing, unexpected
and engineering
o
behavior in a simulation over time. SD is nowadays used as a tool mainly in socio-econom
applications. For an introduction to SD and its applications we refer to Sterman [1].
Simulation runs with varying input parameters can be used to test and evaluate different policies.
absct of meaningful input parameters from the model is selected. They can either be constant
parameters which do not change over time, or they can also be time-dependent control functions. When
stem, it is natural to ask, which of the
two solutions is better? To answer this question properly, one needs to introduce an objective or goal
function. This function is a mapping of the simulated trajectories to a single real number. In this way,
two different simulation runs become comparable: The higher the objective function value, the better. It
is natural to ask, what are the parameters that lead to the best-possible objective function value? This
question cannot be answered from just testing and simulating any finite number of different parameters,
no matter how large this number would be.
Here a
two different runs lead to a different dynamical behavior of the
A System Dynamics model together with control functions and a real-valued objective function is called
a System Dynamics Optimization (SDO) problem in the sequel. The need for an integration of optimization
in the past. In previous approaches, different
methods are used, such as nonlinear local optimization (for example, gradient search methods) or heuris
(such as genetic algorithms), which are essentially based on an “optimization through repeated simulation:
approach, see [2-9]. All these approaches have in common that they at best only provide feasible or local
optimal solutions. Moreover, as pointed out by Burns and Janamanchi [10], nonlinear optimization methods
rely on the availability of derivative information from sufficiently smooth functions. This restriction is in
conflict with certain SD modeling elements, such as if-then-else clauses. In principle, all these methods
methods into the SD framework has been recognized alread;
cannot give a proof of global optimality or any other quality certificate of the computed solutions. The
ultimate goal of our research is to fill this gap and develop a computational method that is able to handle
all kinds of SD models and does not only yield feasible solutions, but also a certificate of optimality.
We formulate the SDO problem as mixed-integer nonlinear program (MINLP). Solving optimization
problems from thi s is theoretically intractable and also known to be computationally difficult in
general. By “solving” we mean to compute a feasible solution for a given instance of the problem together
with a computational proof of its optimality. Therefor we apply the general framework of a branch-and-
bound approach, where the bounds are obtained from relaxations of the original model. To this end, we
relax the MINLP first to a mixed-integer linear program (MILP) and then further to a linear program (LP),
which is solved efficiently using Dantzig’s simplex algorithm [11]. The so obtained solution value defines
a (lower) bound on the optimal value of the original MINLP problem. In case this solution is MINLP
feasible, it would be a proven global optimal MINLP solution. However, this rarely happens in practice.
Hence we either add cutting planes to strengthen the relaxation, or we decide to branch on a variable. As
an example, consider the nonlinear constraint f : 2 +> 2|2|In the LP relaxation this function is replaced
by a polyhedral (linear) outer approximation, which is iteratively refined during the branch-and-bound
process by branching on variables (spatial branching), see Figure 1. For more details on cutting planes
and branch-and-bound for MILP we refer to Nemhauser and Wolsey [12], and for an application of this
framework to global mixed-integer nonlinear programming to Smith and Pantelides [13], and Tawarmalani
and Sahinidis [14,15]. Information on the MINLP framework SCIP which we apply is given in Achterberg
[16], and in particular on nonlinear aspects of SCIP in Berthold, Heinz, and Vigerske [17].
A
i
a) b)
Figure 1: a) Polyhedral outer approximation of f : 24 x
spatial branching.
|, b) initial spatial branching on zero, c) further
As a modeling language, mixed-integer programming offers a high flexibility in formulating various
types of SD modeli As an example, we demonstrate this flexibility by reformulating an if
Say, the SD models contains the following constraint: if « > a then y = 1 else y = 0
me that there are lower and
then-else clause
end. Here x,y are variables to be determined in the simulation. We a
upper bounds known for «, that is, x € [x,Z]. The (binary) variable y is activated (y = 1), if the stock x
becomes larger than a certain given quantity a. Otherwise, y is inactive (y = 0). Note that this condition
constitutes a discontinuous relation between x and y. Due to this discontimuity, an optimization method
based on derivatives would not be applicable (or, to be more precise, would only be applicable after some
smoothing tricks, that come along with further numerical issues). In a mixed-integer formulation of this
condition we introduce two auxiliary variables €+,€- > 0, and then state the following constraints:
g-a = €t-£&, (1a)
F< (©-a)y, (1b)
€& < (a-z)1-y), (Ie)
y € {0,1}. (1d)
Note that constraints (1a), (1b), (Ie) are linear (in-) equalities. Hence Dantzig’s simplex algorithm [11]
is able to deal with them. The integrality constraint (1d) is then relaxed to its continuous counterpart
0 <y <1. If in an optimal solution to this LP relaxation the variable y has a fractional value, one can
the integrality condition either by a cutting plane approach or by branching on y. For details
we refer to the literature [12, 18].
re-introdu
on these two approache
2 Our Approach
In the following we describe our approach to solve mixed-integer dynamic optimization problems that
originate from System Dynamics Optimization (SDO) models. As a foundation we use the linear and
nonlinear programming based mixed-integer nonlinear solver SCIP [16]. This framework already handles
the input of an instance, the set up of the model’s variables and constraints, and the overall branch-and-cut
solution process. Initially, we tried to solve SDO model instances using SCIP as a black-box solver off the
shelf. Although SCIP is a highly capable general-purpose solver, the results in our case were devastating.
It became clear that we had to enhance SCIP to detect and exploit the special structure that is inherent to
SDO problems. That means, the time discretization leaves its traces in the variables and the constraints,
and it pays off to use the natural (time) order of the variables to speed up the computation process.
Certain general-purpose preprocessing techniques also needed to be specially tailored towards SDO. These
issues are described in detail in Section 2.2. With their help it is possible to substantially reduce the upper
bound of the linear relaxation (in case of a maximization problem which we consider here without loss of
generality).
2.1 Problem Formulation
We formulate an SDO problem as a dynamic MINLP problem of the following form:
2Ys2)s (2a)
st. 0= f(ajs1,2j,yj,2)), 7 €{0,1,...,T-1}, (2b)
0 < g(xj,4j,25), 7 €{0,1,...,T}, (2c)
2, <a; < j €{0,1,...,T}, (2d)
¥, Su S9j, F€{0,1,...,T}, (2e)
2,54 <3, j€{0,1,...,T}, (2£)
ER, FE (Ol ansgT}; (2g)
yj €Z!x RP4, je {0,1,...,T}, (2h)
2,)€Z°xR"™*, 7 €{0,1,...,T}. (2i)
Here x = (j);=0,1,...7 stem’s dynamic variables for all time steps j. Further, y
denotes the static variables at each time step j, and z = (z;)j=0,1....,.7 are the control variab!
step j. A variable is said to be dynamic, if its time derivative occurs in the model. A static variable is an
auxiliary variable that aggregates informations from other variables subject to a functional relationship. A
ive function c assigns a real value
,2), and thus allows to compare two solutions with each other. The goal i
identify a global optimal solution with respect to the measure of this function c. Without loss of generality
control variable models the external influence to the system. The obj
to each feasible solution
we focus on maximization problems in equation (2a). The system’s dynamical structure is described by u
continuous nonlinear functions f : R" x R" x Z4 x RP-4 x Z* x R"-* — R". Equations (2b) couple the
dynamic variables of time step j +1 with all (static, dynamic, and control) variables of time step j. These
equality constraints are used to describe discretized differential equations. The v continuous nonlinear
functions g : R" x Z4 x RP-4 x Z* x R’-* — RY in equation (2c) model static relations for a single time
step j. Similar to the distinction of static and dynamic variables, we say that a constraint is dynami
if
it contains (discretized) derivates, and static otherwise. Static constraints can be inequality constraints
ity constraint is covered by using two inequalities of opposite
direction). Bounds on the range of the state variables are given for all time steps by constraint (2d)
in general. (The special case of an equ
for the dynamic variables and by constraints (2e) for the static variables. The dynamic variables are all
continuous variables as specified in equation (2g). Some of the static and control variables may carry
additional integrality constraints, which is specified in equation (2h) and equation (2i), respectively.
2.2 Bound Propagation
Clearly, the bound constraints (2d), (2e), and (2f) are special cases of the static constraints (2c). Never-
we stated them explicitly in our problem formulation (2), bec
se of our solution approach. The bound values 2,2 that are given in (2) are usually weak in
the following sense: They reflect the modeler’s idea on the domain that the optimized system should never
leave. It is known that already “normal” SD models (i.e., without the aspects of control and optimization)
can have a surprising, unexpected behavior due to the nonlinear structure and feedback loops. This is even
more true for SDO problems
tight bounds on the variables is practically impossible under these circumstances. So it is more likely to
expect that initially some rough estimation will be given.
A crucial step in solving mixed-integer linear or nonlinear problems is the preprocessing phase. In this
initial step of the solution process, one tries to extract as much information as possible from the problem at
hand. These information are implicitly contained in the model formulation, but hidden for the solver. The
goal is to make them visible, so that the solver can use them to shorten the solution process. Preprocessing
or presolving refers not just to a single method, but a whole bunch of various techniques, where some are
more general, and others are more specific for certain problem struc
surveys to MILP and MINLP solution techniques including preproce
Bound propagation belongs to the class of general techniques
linear constraint of the form )7/_, a;w; < b, where a;,b are constant
constraints w, < w; < 7; for i= 1,2
ase they are crucial in the initial
where the unintuitive dynamic interacts with an external control. Giving
ures. For more details and general
ing, we refer to [18-22].
ume that the problem contains a
and w; are variables, and the bound
;n. Then for any j with a; > 0 we can try to update the upper
bound via
b aj aj
Wj,—— 7 wi Wi >, (3)
a a
9 ixja>o % iAj,ai<0
and for any j with a; < 0 we can try to update the lower bound via
a ay
- > —w;— > “Dip. (4)
i#jai>0 By if#j.ai<0 od
of a continuous nonlincar constraint f(w1,..., wn) <b with f IR" > R, and bounds w, < wi < 7%
1,2,...,n one needs to solve an unconstrained optimization problem to obtain improved bounds
min {w;,max{f(wy,..-,wp) | w; < wi < W;,%
wn) |W; < wi STF = 1,2,..
(a)
(6b)
max {w,,min{ f(wi,
The computation steps in equations (3), (4), and (5) are repeated for all variables j and all constraints of
the model, until no further bound improvement is achieved.
Bound strengthening is an important step for two reasons. First, we apply linear programming to
solve the subproblems in the branch-and-bound tree. If the variable bounds are weak, then the solution
of the linear programs will also yield weak (upper) bounds for the solution of the MINLP. Second, when
branching on a variable, it needs more branching decisions to fix a variable when the bounds are weak.
Both issues are related to the solution speed and the memory requirement during the solution process.
Hence in general one is interested in obtaining good (narrow) bounds on all variables, assuming that the
time spend for computing them is marginal compared to the gain in solving the overall problem.
In Fiigenschuh et al. [23] it was noted that this gencral bound strengthening procedure gives weak results
on MILPs with a dynamic structure. In this work, a finite difference approach for a partial differential
(transport) equation on a network was examined. It was concluded that the dynamic structure must
be exploited in order to speed up the presolution process by some orders of magnitudes. Furthermore,
dealing with a single constraint at a time yields weaker results compared to the case of multiple constraint
preprocessing. The situation for dynamic MINLP problems approximating ordinary differential equations
s somehow similar. Again, we use the induced structure of the time axis for an ordering of the preprocessing
steps. And then, we do not use a single constraint in each iteration, but certain larger portions of the
model in order to arrive at better (more narrow) bounds. The detail
below.
We define the problem P(w;h,k) for a variable w € {01,...,2n,Yts---sYprZ1s.0-s2r}s and hy ke =
0,1,...,7 with h—k>0as
of this procedure are described
P(z:h,k) = max w, (6a)
st. 0= f(tj41,2j,4j,%), FE{h—k, (6b)
0 < g(xj,yj,2)), JE {h—k,. (6c)
2; $a; <7, je {h—k,... (6d)
y, Su S je f{h—k, (6e)
je{h—k (6f)
je {h—k,...,h}, (6g)
yy €Z2 x RP-4, je {h—k,...,h}, (6h)
2,€Z°xR"*, je{h—k,...,h}. (6i)
Similar, we define the problem P(w;h,k) := min w, with the same constraints as before in (6). These
values are used to update the upper and lower bounds on variable w, respectively. These two auxiliary
problems are very similar to the original problem (2), where the main differences are the exchange of the
objective function, and the smaller time horizon, from h — k to h instead of 0 to T.
We use the solution of (6) in the following way. We select a value for the look-back presolve level k.
For smaller values of k the solution of (6) is faster, but also the obtained bounds are weaker, compared
to large values of k. Here a trade-off between solution speed and a quality of the results needs to be
established. For each h = k,k +1,...,T and each variable z € {21,..-,2n,Yts-++sYps21s++++2r} we solve
for P(w;h,k) and P(w;h, k), where these values are compared with the previously existing best bounds
for w, and updated if necessary. These runs are in principle independent of each other. However, it pays
off to solve in upwind direction of the time, that is, first for h = k, then h = k +1, and so on. The
bounds obtained in earlier runs can directly be re-used in subsequent iterations, and speed up the solution
. Moreover, after one single run all bounds already obtained their respective final values, and it is
ary to repeat this run a second time (this would not lead to better bounds).
3 Test Instances
As test instances we use two system dynamic models that we describe in the following. The models are
taken from the literature (c.f. Bossel [24-26]). We add a control function to each model that allows for an
external interception. The question is how to chose this control function with respect to some optimality
criterion. A typical goal is to interfere with the models’ inherent dynamics in the least possible
then transform each model into a nonlinear mixed-integer problem. In this form it can be used as input
data for our solution algorithm SCIP and our presolve methods.
We
3.1 Predator-Prey
As a first test instance we use a predator-prey system with capacities. The math
formulation of predator-prey models can be dated back to the work of Lotka [27] and Volterra [28]. The
basis of our model’s description below are taken from Bossel [25]. Since Bossel describes a simulation
problem only, we add a control function to the model: We aim to maximize the amount of prey that
can be removed from the system without inducing its collapse. This can be thought of as a model for
sustainable human hunting activity. Since we have to deal with a finite time horizon, it needs to be assured
that not all prey is killed at the (somehow artificial) end-of-time, in order to simply maximize the objective
function.
The SDO diagram is shown in Figure 2. In the mathematical model cor g to this figure, we
introduce three time-dependent functions: the prey population x(t), the predator population y(t), and the
ee — edt (Y
Ye
Oe SK
—_—_ A
D>
ystem Dynamics diagram for the predator-prey control problem.
Figure 2
control z(t). If both populations, predator « and prey y, would be independent of each other, the prey
population would grow proportional to its size and net growth rate (because of the assumed unlimited
capacity). This would lead to an exponential growth with rate a, ie., “ = ax. Without prey, the
predator population would decrease with an exponential rate —d, i.c., 4 ay When both populations
share the same habitat, they meet cach other, and as a result, parts of the prey are consumed by the
predators. The chance of a meeting is proportional to both population sizes xy. In case of a meeting the
prey population is reduced by a factor b and the predator population increases by a factor c. The prey
is additionally hunted by humans. We assume that these hunters can take out any desired part of the
population, independent from its current size. ‘Thus it is reduced by = in each time step. Hence we have
the differential equation s ax — bry — z and % y — dy. As boundary values, the initial
population for predator and prey is given: (0) . Furthermore, the prey population is not
allowed to fall below a certain limit over the whole time horizon: x > €. This in particular prevents to
take out all remaining prey in the very last time step of the model.
We apply a forward discretization to the differentials. Let (At) denote a time discretization. We
introduce continuous variable 2;, yi, 2; € R for i = 0,1,...,7. Then the discretized SDO problem is the
following:
(7a)
= ax; —beiy; —%, 1=0,1,. (7b)
a = -dy; tery, i=0,1,. (7c)
>é i=0,1, (7d)
29 = 27, Yo = Yr, (7e)
= T,Yyo = 4, (7£)
i, Yi % ER, oT. (78)
We emphasize equation (Ze). These constraints enforce any feasible solution to the problem being periodic,
that is, the final state must equal the initial state of the system. Then it is guaranteed that a sustainable
(optimal) solution is reached, which can be continued forever.
3.2 Mini-World
The second test instance is a mini-world model, introduced by Bossel [26]. It is a simplification of the
much more sophisticated World? and World3 model of Forrester [29] and Meadows et al. [30]. World3
uses 18 stocks, 60 parameters, 52 tabular functions, and around 200 equations to model the system’s
relations. Bossel’s Miniworld is an aggregated version, that comes with just 3 stocks: the world population,
the production (industrial, commercial, agricultural), which equals the consumption of goods, and the
environmental pollution. Ca ly, the number of parameters, tables, and equation relations are
much lower. Interestingly, the model shows a qualitatively similar behavior compared to the much more
evolved World3 model. Again, Miniworld is designed as a pure SD model for simulation. We add a control
function to obtain an SDO problem. Our goal is to maximize the economic growth level, in order to provide
a high standard of living to the world’s population. However, if this standard would be too high, then a fast
growing population would quickly consume its natural resources, and thus the population would collapse
soon after. To prevent such behavior we introduce a population level € as a lower limit. The question we
want to address is, how much sustainable growth can bear this mini-world at most for a population being
at least €.
The SDO diagram is shown in Figure 3. The corresponding mathematical model has four time-
dependent functions: p(t) for the population, e(¢) for the envi al pollution, ¢(¢) for the ¢¢ i
and 2(t) for the control function (the growth level).
We apply a forward discretization to the differential terms, with (At) being the size of a step in the
time discretization. The continuous variables p;,¢;, ci, 21 € R for i = 0,1,..., 7° approximate the functions
of the respective same name. Furthermore, we introduce the following real-valued continuous variables,
each for i =0,1,...,T: a; are the number of births and Q; are the number of deaths at time step i. 7% is
the consumption level. 1; des
time, and for this we introd:
bes the environmental quality. The environmental conditions change over
se €; for the environmental recovery and ¢; for the environmental destruction.
Using a forward Euler schema, the discretized SDO problem reads as follows:
(8a)
PF Ly 1, FHI oT (8b)
a; = 0.03-pi-mi-yis 1=0,1,...,7, (8c)
Q; = 0.01-p;-e:, 1=0,1,...,T, (8d)
Se seis i=0,1,...,T7, (8e)
G =0.02-p;-71, 1=0,1,...,T, (8f)
e,=ef —e7 +10, i=0,1, (8g)
ef BQO. TH Os lywany (8h)
e7 <10-a, i=0,1,...,T, (8i)
e=01-(10-—e7), #=0,1,...,7, (8))
10=e-:, 1=0,1,...,T, (8k)
Ci41 — Ci
a (81)
a=%, (8m)
p2g, t=tt+1, (8n)
Po =P, (80)
eo =6, (8p)
oo =%, (8a)
Di €i,e} .e; (8r)
x; € {0,1}, (8s)
The change in the population p from time step i to i+ 1 is a result of the births a; minus the deaths
Q, in time step i (8b). The number of births a; in time step i is proportional to size of the population
pi, the environmental quality jz;, and the consumption level 7;, where 0.03 is the proportionality factor
(birth rate) (8c). The number of deaths ; in time step i is proportional to the population p; and the
environmental pollution ¢;, with a proportionality factor (death rate) of 0.01 (8d).
> Env. quality (u)
regeneration (e) Pa Destruction (2)
raion) | —§<—<—>
i \ {ft
\\
< Ne
NI
Figure
em Dynamics diagram for the mini-world control problem.
The change in the environmental pollution e from time step i to i-+1 is a result of the negative environ-
mental destruction ¢; and the positive environmental recovery ¢; in time step i (8c). The environmental
destruction ¢; in time step i is proportional to the population size p; and the consumption level 7; (8f).
Here 0.02 is the proportionality factor (environmental destruction rate). The environment is able to recover
over time. If the environmental quality j1; is above a certain threshold value (here 1.0), then the recovery
e; in time step i is proportional to the environmental pollution ¢;, with 0.1 being the proportionality factor
(recovery rate). However, if the environmental quality j; is below the threshold value, then the recovery
rate is no larger than 0.1 (8j). For a proper representation of this constraint in our model, we introduce a
binary variable x; for the decision if ¢; is above or below the threshold value of 1.0 (8g). The part above
and below is separately stored in two auxiliary variables e* (8h) and e7 (8i), respectively. Note that the
variable jx; docs not appear in the equations (8j) directly. Equation (8k) reflects the dependency of the
environmental pollution and its inverse, the environmental quality. This equation is responsible that 4; is
indirectly part of equations (8}).
The change in the production and consumption c from time step i to i+1 is the result of a logistic growth
function of Verhulst type (81), which depends on both the consumption level ; and the environmental
pollution e;. The control variable 2 plays the role of the system’s capacity (this is a constant in the original
Verhulst equation). The constant value of 0.05 is a growth rate for the consumption. The consumption
level 7; equals the production ¢; (8m).
The population p; must not fall short of the given level &; in each time step i (8n). Initially, in time
step i = 0, the size of the population (80), the environmental pollution (8p), and the consumption resp.
production (8q) are given.
4 Computational Results
We discuss the application of our methods to the two test problems. We start with the smaller and thus
simpler predator-prey problem, and present the mini-world results afterwards.
4.1 Predator-Prey
Figure 4 shows the population sizes for different levels of constant hunting activity. That is, 2 for
all i = 0,1,...,7, and 2 € [0,0.7]. For 7 = 0, ic., no hunting activity, this reduces to the classical
Lotka-Volterra model. The predator and prey populations (dark blue lines) follow a cyclic behavior with
a period length of 7 time units (say, years). If positive hunting activity comes into play, i.c., 7 > 0, then
the prey population will sooner or later become extinct. This can be seen in detail in Figure 5, where even
nall positive but constant hunting activity leads to an extinction after a few We conclude that
a constant hunting level is not sustainable for this dynamical system.
a
Predator
bp oR
Population
Population
Time
Figure 4: Trajectories for the predator and the prey for different levels of a constant control, varying in
(0, 0.7].
level ScIP_0 10 20 40 80 120
time [min] [0.005 6 +9 13 52 144 200
Table 1: Running times for different presolve look-back levels for the predator-prey SDO.
the bounds on the variables need to be stength-
subroutine determines the time spent here, as
In order to apply mixed-integer nonlinear programmin,
ened in a preprocessing step. The look-back level of this
well as the quality of the result. In general, a higher look-back level leads to higher running times, but
also better results. In Figure 6 we compare the remaining bounds of the original presolve routines of the
MINLP solver SCIP with our results for various look-back levels. In Table 1 we show the corresponding
running times.
Figure 7 shows a solution to the predator-prey problem with human hunting activity. Using our presolve
methods we could prove that this solution is at most 0.77% away from the lower bound. This means, this
solution might be optimal, but in case it is not optimal, an optimal solution can only be 0.77% better.
We have a high level of hunting for about 4 years, and a close season for about 3 years allowing the prey
population to regenerate. At the end of the 7 years cycle both populations are of the same size as in the
beginning. This means that the same trajectories can be repeated to infinity, hence a sustainable control
schema of the s
In Figure 8 we compare three phase diagrams in the predator-prey space. In black we show the
population dynamic without hunting activity. The blue curve is the population dynamic for constant
hunting activity. The red curve is the population dynamic for the sustainable hunting schema.
tem is achieved.
4.2 Mini-World
We further simplify the control z. To this end, we introduce a piecewise constant control with two states
zo and 21, where zo is valid for the first T/2 many time steps, an 2 is valid for the second half of T/2
time steps. Then the objective function is a two-dimensional function depending on the values z, 21 for
the two constant controls. Figure 9 shows a contour plot of the objective function values for the two
controls varying independently in the interval [1,10]. Feasible solutions are in the lower-left corner. The
Population
15
Time [years]
Figure 5: Trajectories for the predator and the prey for constant hunting level of 0.01. After a few cycles,
the prey population becomes extinct.
—— Standard bound propagali
1g -|| $0 bound propagation,
—— $D bound propagation,
Predator Population
8 oe
=
. sf spe
a.
4
2
0 1 6 7
3 4
Time [years] ( t=0.01)
Figure 6: Presolved bounds for different look-back levels.
black curve marks the borderline between the feasible and the infeasible region. Our goal is to keep the
mini-world population above € > 4 (billion inhabitants). Any solution having less inhabitants is infeasible.
Solutions in the infeasible region all fall short of the population lower limit of constraint (8n). A possibly
optimal solution occurs for 29 = 1.79 and 2, = 2.16. This solution yields the highest total consumption
level, while keeping the population size above the specified limit.
In Figure 10 we show three different solutions. The first for 2) = 21 = 1 is feasible but not optimal.
The second for zo = 21 = 10 is infeasible. (These two reference plots can also be found in Bossel [26].)
Our potential optimal solution for z) = 1.79 and 2 = 2.16 is shown as third plot.
Note that the plot in Figure 9 and the identification of an optimal solution is based on a number
of simulation runs. Even if we use a fine grid, we cannot be absolutely sure that there is not a better
solution that is hidden between two neighboring grid points. Our goal is to demonstrate that by using
our bounding techniques we can cut off large portions of the search space. It is thus guaranteed that no
better solution can be found in such a cut off part. We apply our bounds strengthening presolve routine.
The results for different look-back levels are shown in Figure 11 and Table 2. We vary the look-back level
between 0 and 10. The row “SCIP” shows as reference values the results of the solver SCIP without using
our bounds strengthening routines. For cach run, we give the CPU running time and the dual bound for
the presolve phase. After the presolve phase is finished, SCIP starts the branch-and-bound phase with a
time limit of one hour. We give the best dual bound and the best feasible solution found in the end, and
10
Population / Control
time
0 03 i 4
Prey
Figure 8: Phase diagrams for predator-prey population without hunting (black), constant hunting (blue),
and sustainable hunting (red).
the gap, g = 42 . 100%, where p is the primal and d is the dual bound. The SCIP-level (first line in
Table 2) is special: The solver was very fast, but due to its numerical instability “proved” that the problem
is infeasible (i.c., does not have a feasible solution), which is of course not true.
A part of the search space that does not require further attention due to the bounding argument is
shown in Figure 12. Here a look-back presolve level of 5 time steps was applied. Our optimal solution for
z = 1.79 and z, = 2.16 has an objective function value of 230.14. Any solution in the box defined by the
green lines has an upper bound on the objective function value of less-or-equal 229.94, hence no better
solution can be found in this part of the control space.
age of our implementation, we are, however, not able to prove that our potential
is indeed a global optimal one. The solver SCIP terminates, even in the highest presolve
level, with a huge gap of 61.26%. Our ongoing research is devoted to close this gap to 0% within the given
time limit, and to solve the problem for a control function with more than two constant segments.
5 Summary and Conclusions
We presented System Dynamics Optimization (SDO) models as an extension of classical System Dynam-
ics Simulation (SD) models, where one (or several) control mechanism are introduced into the s)
Moreov specified that the modeler wants to achieve. This raises a theoret
challenging question how to actually compute such control function
stem.
_ a desired state must be
well as prac
as
and in case there
Figure 9: Objective function values for piecewis
optimal solution is in (1.79, 2.16).
onstant controls varying in (1, 10] x [1, 10]. The (unique)
Constant Low Control Constant High Control Optimal Control
ee
6 6
Ba Ba
g 8
wn wn
2 2 as
0 0
0 50 100 0 50 100
Time Time
a) b) °)
Figure 10: a) Feasible but suboptimal solution with control z9 = 21 = 1. b) Infeasible solution with control
29 = 21 = 10. c) Optimal solution with zp = 1.79 and z; = 2.16.
is more than one possible control, how to identify the globally best among all possibles, with respect to
some optimality criterion. To answer this question one can simply resort to an optimization-by-simulation
approach, where a certain number of simulation runs is carried out, and the control functions vary between
each run. The mechanism to alter the control functions can be found in the literature, for example, Ge-
netic Algorithms, Tabu Search, or Nelder-Mead Simplex. We demonstrated that these methods all have
the principle disadvantage that they do not give any kind of performance guarantee on the computed
result.
To overcome this weakness, we suggest mixed-integer nonlinear programming (MINLP) methods. To
apply this method, the ordinary differential equations in the SDO problem are approximated by some
discretization schema. Certain features that are common in the SD modeling language (such as tabular
functions, or if-then-else constructs) also need to be transcribed and reformulated by introducing piecewise
linear inequalities and binary or integer variables. Then this MINLP model can be relaxed to a linear
programming problem, and solved by Dantzig’s simplex method. This results in a bound on the objective
function value to the SDO problem, which is then improved by adding suitable cutting planes and a
branch-and-bound procedure. In order to achieve good bounds quickly, we further enhance
using a specially tailored bounds strengthening as preprocessing, and a primal heuristic to identify feasible
solutions earlier.
We applied our method to two test problems from the literature. We showed that our bounds strength-
ening procedure significantly redu h space, which directly relates to the time spent to compute
an optimal solution. In the predator-prey model with additional control by human hunting, we were able
1 = sentna bond poner || if pe
9 |] ——50 bound propagation J f
—— SD bound propagati 4 a
8-| —— 5p bound propagation, eae
a yy
5
Z
a
0 10 20 30 40 50 60
Time [years](_ t=0.50)
Figure 11: Presolved bounds for different look-back levels.
presolve branch-and-bound
level total time [h| dual time [h| dual primal gap
SCIP 0.002 = = = = =
0 1.33 6.253-10° 1 6.253-10° 184.655 > 100%
ap 4.42 6.395-10? 1 6.37910? 184.655 99.4
2 2.702 5.806-10? 1 5.787 +10? 184.655 68.1
5 4.53 5.082-10? 1 5.062-10? 184.655 63.5
10 9.01 4.78610? 1 4.766- 107 184.655 61.3
Table 2: Computational results for different presolve look-back levels.
to compute a control schema that guarantees a maximum yield of prey under the condition that we seck
a sustainable solution, so that neither predator nor prey are dying out at any time in the future. The
optimal control has both hunting season and a close season, whereas any constant schema that comes
without such close season, no matter how little prey is hunted, will lead to an extinction of both species.
For the mini-world control model we demonstrated that larger areas of the solution space do not need
to be examined by a bounding argument, once a suitable good (perhaps optimal) solution is found. We
computed a piecewise constant globally optimal control that achieves a maximum consumption level under
the additional constraint that humankind shall not be perished from this earth.
We believe that our methods are generally applicable to all different kind of SDO models. At the
nt stage of our work, each model needs to be individually analyzed and treated. We are working
pre
towards a black-box solver that a model will be able to use for his or her models in the same way SD
simulation tools are used today. We think that having not only good-looking solutions, but additionally a
certificate of optimality, is a further asset in practical applications.
Acknowledgement
We gratefully acknowledge the funding of this work by the German Research Association (Deutsche
Forschungsgemeinschaft, DFG), Collaborative Research Center CRC1026 (Sonderforschungsbereich SFB1026).
230
26
225
24
220
22
230,14 215
i 210
Fie 2 A
205
1.64 SS 200
1.4} 195
bo
12 w 190
7 SSS 185
1 25
2(0)
Figure 12: The part of the control space with 2; < 1.3 and 29 < 1.3 marked by the gray box can be cut
off. No better feasible solutions can be found in here.
References
1]
Sterman, J. D., 2000. Business Dynamic
Irwing McGraw-Hill, Boston.
~ Systems Thinking and Modeling for a Complex World.
Duggan, J., 2008. Using System Dynamics and Multiple Objective Optimization to Support Policy
Analysis for Complex Systems. In H. Qudrat-Ullah, J. Spector, P. Davidsen, editors, Complex Decision
Making, Understanding Complex Systems. Springer Verlag, Berlin, 59-82.
Elmahdi, A., Malano, H., Etchells, T., Khan, $., 2005. System Dynamics Optimisation Approach
to Irrigation Demand Management. In Proceedings of the MODSIM 2005 International Congress on
Modelling and Simulation. 196-202.
Dangerfield, B., Roberts, C., 1996. An Overview of Strategy and Tactics in System Dynamics Opti-
mization. Journal of the Operational Research Society, 47:405-423.
Kleijnen, J., 1995. Sens
analysis and statistical design of experiments. §
vity analysis and optimization of system dynamics models: Regression
stem Dynamics Review, 11(4):275-288.
Fu, M., 1994. Optimization via simulation: A review. Annals of Operations Research, 53(1):199-247.
Keloharju, R., Wolstenholme, E., 1989. A Case Study in System Dynamics Optimization. Journal of
the Operational Research Society, 40(3):221-230.
Keloharju, R., Wolstenholme, E., 1988. The basic concepts of system dynamics optimization. Systemic
Practice and Action Research, 1(1):65-86.
Gustafson, L., Wiechowski, M., 1986. Coupling DYNAMO and optimization software. System
Dynamics Review, 2(1):62-66.
Bums, J. R., Janamanchi, B., 2007. Optimal Control and Optimization of System Dynamics Models:
Some Experiences i lations. In Proceedings of the 2007 Meeting of the Southwest Region
Decision Scienc:
Rec 2
Institute.
Dantzig, G., 1963. Linear programming and extensions. Princeton University Press and RAND
Corporation.
[12]
[13]
[14]
[15]
[16]
f17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
rx)
x
[26]
27
[23]
[29]
[30]
Nemhauser, G. L., Wolsey, L. A., 1988. Integer and Combinatorial Optimization. John Wiley & Sons,
New York.
Smith, E., Pantelides, C., 1999. A Symbolic Reformulation/Spatial Branch-and-Bound Algorithm for
the Global Optimization of nonconvex MINLPs. Computers & Chemical Engincering, 23:457-478.
Tawarmalani, M., Sahinidis, N., 2004. Global optimization of mixed-integer nonlinear programs: A
theoretical and computational study. Mathematical Programming, 99:563-591.
Tawarmalani, M., Sahinidis, N., 2005. A polyhedral branch-and-cut approach to global optimization.
Mathematical Programming, 103:225-249.
Achterberg, T., 2009. SCIP: Solving constraint integer programs. Mathematical Programming Com-
putation, 1(1):1-41.
Berthold, T., Heinz, $., Vigerske, S., 2012. Extending a CIP Framework to Solve MIQCPs. In J. Lee,
S. Leyffer, editors, Mixed Integer Nonlinear Programming, volume 154, part 6 of The IMA Volumes
in Math and its Appl . Springer Verlag, Berlin, 427-444.
Fiigenschuh, A., Martin, A., 2005. Computational Integer Programming and Cutting Planes. In
R. W. K. Aardal, G. Nemhauser, editor, Handbook on Discrete Optimization, Series Handbooks in
Operations Research and Management Science, volume 12. Elsevier, 69 — 122.
Belotti, P., Lec, J., Liberti, L., Margot, F., Wachter, A.,
techniques for non-convex MINLP. Optimization Metho
2009. Branching and bounds tightening
and Software, 24(4-5):597-634.
Belotti, P., Kirches, C., Leyffer, $., Linderoth, J., Luedtke, J., Mahajan, A., 2012. Mixed-Integer
Nonlinear Optimization. Technical Report ANL/MCS-P3060-1112, Argonne National Laboratory,
Argonne, Illinois.
Burer, $., Letchford, A.
Surveys in Operations Res
2012. Non-convex mixed-integer nonlinear programming: A survey.
rch and Management Science, 17(2):97-106.
Savelsbergh, M., 1994. Preprocessing and Probing Techniques for Mixed Integer Programming Prob-
lems. INFORMS Journal on Computing, 6:445-454.
Dittel, A., Fiigenschuh, A., Géttlich, $., Herty, M., 2009. MIP Presolve Techniques for a PDE-based
Supply Chain Model. Optimization Methods and Software, 24(3):427 — 445.
Bossel, H., 2007. Systems and Models: Complexity, Dynamics, Evolution, Sustainability. Books on
Demand GmbH, Norderstedt.
Bossel, H., 2007. System Zoo 2 Simulation Models: Climate, Ecosystems, Resources. Books on
Demand GmbH, Norderstedt.
Bossel, H., 2007. System Zoo 3 Simulation Models: Economy, Society, Development. Books on
Demand GmbH, Norderstedt.
Lotka, A. J., 1925. Elements of Physical Biology. Williams and Wilkins Company.
Volterra, V., 1926. Variazioni e fluttuazioni del numero d’individui in specie animali conviventi.
Memorie della Regia Accademia Nazionale dei Lincei. Ser. VI, 2:31 — 113.
Forrester, J., 1973. World Dynamics. Waltham, MA, Pegasus Communications.
Meadows, D., Meadows, D., Randers, J., Behrens III, W., 1972. The Limits to Growth. Universe
Books, New York, NY.
