18th IEEE International Conference on Control Applications

Part of 2009 IEEE Multi-conference on Systems and Control

Saint Petersburg, Russia, July 8-10, 2009

A Generic Dynamic Programming Matlab Function

Olle Sundström and Lino Guzzella

Abstract— This paper introduces a generic dynamic pro- The class of optimal control problems that can be solved
gramming function for Matlab. This function solves discrete- using the proposed Matlab function can be written as:
time optimal-control problems using Bellman’s dynamic pro-
gramming algorithm. The function is implemented such that the min J(u(t)) (1)
user only needs to provide the objective function and the model u(t)
equations. The function includes several options for solving s.t.
optimal-control problems. The model equations can include
several state variables and input variables. Furthermore, the ẋ(t) = F (x(t), u(t), t) (2)
model equations can be time-variant and include time-variant x(0) = x0 (3)
x(tf ) ∈ [xf,min , xf,max ]
state and input constraints. The syntax of the function is
explained using two examples. The first is the well-known (4)
Lotka-Volterra fishery problem and the second is a parallel x(t) ∈ X (t) ⊂ R n
(5)
u(t) ∈ U(t) ⊂ Rm
hybrid-electric vehicle optimization problem.
(6)
I. INTRODUCTION
where
When developing causal suboptimal controllers it is an
J(u(t)) = G(x(tf )) + ∫
tf
advantage if the optimal controller is known, even if this con- H(x(t), u(t), t)dt (7)
troller is not causal. In many cases, such optimal controllers 0

can be found using the deterministic dynamic programming is the cost functional. The important characteristics of the
(DP) algorithm introduced in [1]. Of course, this optimal considered optimal-control problems are the time-variant
controller can be found only if all future disturbance and constraints on the input and the state, the constrained final
reference inputs are known. In this sense, this solution is state, and the time-variant model equation. These problems
not causal. Nevertheless, this optimal solution is very useful, are in general difficult to solve. However, these problems
because it can be used as a benchmark to which all other can be solved using Bellman’s DP [1], provided that the
causal controllers can be compared to. conditions stated in the next section are satisfied.
Many excellent text books have been published on the
II. DYNAMIC PROGRAMMING ALGORITHM
subject of DP theory, among them [2] and [3]. An overview
of the history and development of dynamic programming This section gives a brief overview of the deterministic DP
is shown in [4]. Interested readers are referred to these algorithm as it is implemented in the dpm function. Since
references for a detailed discussion of the basic ideas of DP is used here to solve a continuous-time control problem,
DP. When implementing the deterministic DP algorithm on the continuous-time model (2) must be discretized in time
a computer there are many numerical issues that arise that first. Let the discrete-time model be given by
have, so far, not yet received sufficient attention. Also, since xk+1 = Fk (xk , uk ), k = 0, 1, . . . , N − 1 (8)
the computational complexity of every DP algorithm is ex-
ponential in the number of states and inputs, special attention with the state variable xk ∈ Xk and the control signal uk ∈ Uk .
must be given to minimizing the overall computational cost. A. Basic Algorithm
Let π = {µ0 , µ1 , . . . µN −1 } be a control policy. Further let
The implementation of suitable numerical algorithms that
efficiently solve a given DP problem is, therefore, a nontrivial
part of a design process. the discretized cost of (7) using π with the initial state
This paper presents a Matlab function that efficiently x(0) = x0 be
solves deterministic DP problems. The focus lies on the opti- Jπ (x0 ) =gN (xN ) + φN (xN ) . . .
mal control of non-linear, time-variant, constrained, discrete- N −1
time approximations of continuous-time dynamic models. + ∑ hk (xk , µk (xk )) + φk (xk ), (9)
One area where such DP tools have been used successfully is k=0
the energy management problem in hybrid-electric vehicles where gN (xN ) + φN (xN ) is the final cost. The first term
[5], [6]. gN (xN ) represents the final cost in (7). The second term
O. Sundström, Department of Mechanical and Process Engineering ETH is an additional penalty function φN (xN ) that can be used
Zurich, 8092 Zurich, Switzerland and Empa, Swiss Federal Laboratories to enforce a constraint on the final state (4). The function
for Materials Testing and Research, 8600 Dübendorf, Switzerland, sund- hk (xk , µk (xk )) is the cost of applying the control µk (xk )
stroem@[Link]
at xk , according to H(x(t), u(t), t) in (7). The state con-
straints (5) are enforced by the penalty function φk (xk )
L. Guzzella, Department of Mechanical and Process Engineering ETH
Zurich, 8092 Zurich, Switzerland, guzzella@[Link]

978-1-4244-4602-5/09/$25.00 ©2009 IEEE 1625

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 19,2025 at 10:47:09 UTC from IEEE Xplore. Restrictions apply.
 for k = 0, 1, . . . , N − 1. The optimal control policy π o is the In particular the syntax is shown for a simple optimal control
policy that minimizes Jπ problem. Since the problem is simple the entire code is
J o (x0 ) = min Jπ (x0 ),
shown and explained. The dpm function can be downloaded
(10) at [7].
π∈Π

where Π is the set of all admissible policies. When solving discrete-time optimal control
Based on the principle of optimality [1], the DP algorithm problems the dpm function is normally called using
evaluates the optimal cost-to-go1 function Jk (xi ) at every
[res dyn] = dpm(fun,par,grd,prb,options);
node in the discretized state-time space2 by proceeding
backward in time: where fun is the model function handle, par is any user
1) End cost calculation step defined parameter structure that is forwarded to the model,
grd is the grid structure, prb is the problem structure, and
JN (xi ) = gN (xi ) + φN (xi ) (11) options is the option structure. The output of the dpm
2) Intermediate calculation step for k = N − 1 to 0 function are normally two structures representing the DP-
output and the signals from forward simulation of the model
Jk (xi ) = min { hk (xi , uk ) + φk (xi ) . . . using the optimal control input map.
uk ∈Uk
Since the DP algorithm is often time consuming,
+Jk+1 (Fk (xi , uk ))} (12) the dpm function can also be used only for for-
The optimal control is given by the argument that minimizes ward simulation when the DP output structure dyn is
the right-hand side of equation (12) for each xi at time index precalculated. This can be very useful when chang-
k of the discretized state-time space. ing the initial condition or when increasing the start-
The cost-to-go function Jk+1 (x) used in (12) is evaluated ing time N0 of the problem. To call the dpm function
only on discretized points in the state space. Furthermore, when the DP output structure is already calculated use
the output of the model function Fk (xi , uk ) is a continuous
res = dpm(dyn,fun,par,grd,prb,options);
variable in the state space which can be between the nodes
of the state grid. Consequently, the last term in (12), namely All the structures in the code above are further explained in
Jk+1 (Fk (xi , uk )) must be evaluated appropriately. There the remainder of this section.
exist several methods of finding the appropriate cost-to-go
Jk+1 (Fk (xi , uk )) such as using a nearest-neighbor approxi- A. Problem
mation or using more advanced interpolation schemes. In the In the problem structure all necessary parameters that
dpm function introduced in this paper, linear interpolation define the problem are given. The important parameters
of the cost-to-go Jk+1 is used. Since the state and the are the time step Ts of the model description and the
input grids are equally spaced, the computational cost of problem length N. Moreover, in the problem structure an
this interpolation is low compared to the cost induced by the optional cell array can be defined, which contains time-
model evaluations. variant information relevant for the problem description. For
The output of the algorithm (11)–(12) is an optimal control example, if the model explicitly depends on the time the cell
signal map. This map is used to find the optimal control array W{1} would contain a time vector with N elements.
signal during a forward simulation of the model (8), starting The corresponding elements in these time-variant vectors are
from a given initial state x0 , to generate the optimal state forwarded to the model function throughout the problem.
trajectory. In the optimal control signal map the control The problem structure can also contain a starting time index
signal is only given for the discrete points in the state space where the forward simulation starts. This can be helpful
grid. The control signal, therefore, must be interpolated when when searching for a time optimal solution. An overview
the actual state does not coincide with the points in the of the problem structure is shown in Table I.
state grid. In general, the complexity of the DP algorithm
TABLE I
is exponential in the number of state and input variables.
PROBLEM - STRUCTURE ( PRB )

III. DPM-FUNCTION
Ts time step (is passed to the model function)
The dpm function solves the discretized version of the N number of time steps in problem (integer that defines
optimal control problem (1)–(7) using the dynamic pro- the problem length)
gramming algorithm introduced in Section II-A. This section N0 (optional) start time index (only used in forward
shows the syntax and commands for solving such problems. simulation)
W{.} (optional) vectors with length N containing time-
1
The terms cost-to-go and optimal cost-to-go are used equivalently variant data for the model
throughout this paper referring to optimal cost-to-go. It is important to note
that the term optimal is used in the sense of optimality achievable under
the numeric errors. B. State/Input Grids and Constraints
2
The following notation is used: xik denotes the state variable x in the
discretized state-time space at the node with time-index k and state-index The grid structure grd contains all the information about
i. xk denotes a (state-)continuous state-variable at time k. the state and input grids and constraints. An overview of the

1626
Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 19,2025 at 10:47:09 UTC from IEEE Xplore. Restrictions apply.
 TABLE III
grd structure is shown in Table II. The grd structure is com-
OPTIONS - STRUCTURE ( OPTIONS )
posed by cell arrays, where there is a cell for each state vari-
able and each input variable. For example, for a problem with HideWaitbar hide waitbars (0/1)
two state variables the grd structure contains grd.X{1}, Warnings show warnings (0/1)
SaveMap save cost-to-go map (0/1)
grd.X{2}, [Link]{1}.lo, [Link]{2}.lo, and so UseLine use boundary line method (0/1)
on. The input grid is used in a similar way depending on FixedGrid (used if UseLine=1) using the original grid as
the number of input variables of the problem. specified in grd or adjust the grid to the boundary
lines (0/1)
Iter (used if UseLine=1) maximum number of itera-
TABLE II
tions when inverting model
GRID - STRUCTURE ( GRD ) Tol (used if UseLine=1) minimum tolerance when in-
verting model
Nx{.} number of grid points in state grid InfCost a large cost for infeasible states (I=1)
Xn{.}.lo lower limits for each state (vector for time-variant or
scalar for fixed) Minimize (optional) minimizing (or maximizing) cost function
Xn{.}.hi upper limits for each state (vector for time-variant or (0/1) default is minimizing
scalar for fixed) InputType (optional) string with the same number of characters
XN{.}.lo final state lower constraints as number of inputs. Contains the character ’c’ if
XN{.}.hi final state upper constraints input is continuous or ’d’ if discrete (default is all
X0{.} inital value (only used in forward sim) continuous).
gN{1} (optional) Cost matrix at the final time (must
Nu{.} number of grid points in input grid be of size([Link]{1}) = [[Link]{1}
Un{.}.lo (optional) upper limits for each input (vector for time- [Link]{2} ... [Link]{.}])
variant or scalar for fixed)
Un{.}.hi (optional) upper limits for each input (vector for time-
variant or scalar for fixed)
the forward simulation of the model when applying the opti-
mal control input map. The dyn structure is associated with
the dynamic programming algorithm, the optimal cost-to-go,
C. Options
and the optimal control input map. When the boundary line
The DP approach can be used for many different prob- method is used the dyn structure also contains the boundary
lem settings and the options structure defines how to use lines (with the states, inputs, and costs). An overview of the
the algorithm. An overview of the options that can be two structures are shown in Tables IV and V.
specified in the options structure is shown in Table III.
TABLE IV
The HideWaitbar options decides if waitbars are shown
DP OUTPUT- STRUCTURE ( DYN )
or not when running the DP algorithm. The SaveMap
option determines if the optimal cost-to-go is saved and [Link] Xo,Uo{.},Jo contains the cost, input, and state for
returned. Note that the memory requirements increase when the upper boundary line
[Link] Xo,Uo{.},Jo contains the cost, input, and state for
SaveMap=1. the lower boundary line
An important option is the UseLine option, which de- Jo{.,.} optimal cost-to-go (indexed by input number and
cides if the boundary line method, introduced in [8], is used time index)
Uo{.,.} optimal control input (indexed by input number and
or not. The boundary line method is very useful for increas- time index)
ing the accuracy of problems with final state constraints. For
more information about the boundary line method readers are
referred to [8]. Note in the actual version of the dpm function TABLE V
it can only be used when there is only one state variable. If RESULTS - STRUCTURE ( RES )

the boundary line method is used, i.e., if UseLine=1, there

X{.} state trajectories
are three additional options Iter, Tol, and FixedGrid C{.} cost trajectory
that must be defined. The options Iter and Tol determines I infeasible vector (problem is not solved if nonzero
the stopping criteria when numerically inverting the model elements)
function. The option FixedGrid decides whether to adjust signals structure containing all the signals that were saved in
the grid to the boundary lines or fix the grid to the definition the model function
in grd.
Finally, the InfCost is the cost of infeasible states
and inputs of the model. When not using the boundary E. Model
line method InfCost is also used to enforce the final The equations describing the model must be implemented
state constraints in (9), with φN (xN ) =InfCost when in a correct format in order to be used with the dpm
xN ∉ [xf,min , xf,max ]. function. To generate a sample function use the command

D. Output dpm('sample_model',Nx,Nu);

The outputs of the dpm function are two structures, namely This command will save an m-function as
res and dyn. The res structure contains the results from sample_model.m with a random model of Nx state

1627
Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 19,2025 at 10:47:09 UTC from IEEE Xplore. Restrictions apply.
 variables and Nu input variables, suitable for usage with In the first example of the Lotka-Volterra fishery problem
the dpm function, which can be used as a template when the entire code necessary to use the dpm function is shown.
developing a new problem description. In the hybrid-electric vehicle example, however, the model
In general the model function should have the format: function contains far too many lines to be included in this
paper. Interested readers can download the complete model
function [X, C, I, signals] = mymodel(inp,par) equations at [7].
where the model input structure inp is generated by the A. Lotka-Volterra Fishery
dpm-function and contains the elements in Table VI. The In order to evaluate the optimal solution by means of DP
structure par can contain any user defined parameters a discrete-time approximation of the continuous-time Lotka-
necessary in the model function. It is important that the Volterra model is used. Using an Euler forward approxima-
model function preserves the size of the inputs to the outputs. tion with a time step Ts = 0.2 days, the discrete-time model
Consequently, the elements inp.X{.}, inp.U{.} and the is
outputs X{.}, C{.}, and I must have the same size. The
structure signals can contain any user defined internal xk+1 = f (xk , uk ) + xk , k = 0, 1, . . . , N − 1 (13)
signals in model. These signals are stored during the forward
where
simulation and returned in the res structure when calling
x2
f (xk , uk ) = Ts ⋅ ( ⋅ (xk − k ) − uk ) .
the dpm-function, Table V. 2
(14)
100 1000
TABLE VI
INPUT- STRUCTURE ( INP ) The state xk ∈ [0, 1000] is the amount of fishes in a lake, the
control signal uk ∈ [0, 10] is the constant fishing rate during
X{.} current states (n+m dimensional matrix form de- one time step. The discrete-time optimal control problem of
pending on the number of inputs and state variables)
U{.} current inputs (n+m dimensional matrix form de- maximizing the amount of fishes caught over a fixed period
pending on the number of inputs and state variables) of time can be formulated as follows:
W{.} current time-variant data (scalar) N −1
∑ −uk ⋅ Ts
Ts time step
min (15)
uk ∈[0, 10] k=0

s.t.
xk+1 = f (xk , uk ) + xk
TABLE VII
MODEL OUTPUTS
(16)
x0 = 250 (17)
xN ≥ 750 (= xf,min )
X{.} resulting states after applying inp.U{.} at
inp.X{.} (same size as inp.X{.}) (18)
C{.} resulting cost of applying inp.U{.} at inp.X{.} xk ∈ [0, 1000] (19)
(same size as inp.X{.})
I set with infeasible combinations (feasible=0, infeasi- 200
ble=1) (same size as inp.X{.}) N= + 1. (20)
Ts
signals structure with user defined signals (same size as
inp.X{.}) To solve this optimal control problem the model function
(13) is implemented in Matlab as:

IV. EXAMPLES function [X, C, I, signals] = fishery(inp,par)

To illustrate the usefulness of the dpm function, two % state update

func = (0.02.*(inp.X{1}-inp.X{1}.ˆ2/1000)-inp.U{1});
examples are discussed below. First, the well-known Lotka- X{1} = [Link].*func + inp.X{1};
Volterra fishery problem [9] is explained and solved using
the dpm function. Of course, there exist an analytic solution % cost
C{1} = -[Link].*inp.U{1};
to the continuous-time Lotka-Volterra fishery problem, and
it is therefore not necessary to use a DP algorithm to solve % infeasibility
it. However, since this problem is simple and is similar to I = 0;
the problems normally solved with DP, it is used as an signals.U{1} = inp.U{1};
example to illustrate the syntax of the dpm function. Second,
an example of an optimal energy management problem for
a parallel hybrid-electric vehicle is solved using the dpm Since the state and input spaces have to be discretized,
function. This problem is well suited for the DP algorithm. the dpm function includes a simple way to define such
Not surprisingly, DP has been used extensively proposed in grids. Let the state variable be limited between 0 and 1000
the literature to solve such energy management problems, and let it be discretized using a step of 5 such that xk ∈
both for comparison to causal controllers and for evaluation {0, 5, 10, ..., 995, 1000}. Also, let the control input variable
of different system configurations. Some examples are [10], be limited between 0 and 10 and let it be discretized with
[11], [12], and [13]. a step of 0.5 such that uk ∈ {0, 0.5, 1, ..., 9.5, 10}. The

1628
Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 19,2025 at 10:47:09 UTC from IEEE Xplore. Restrictions apply.
 for an initial state of x0 = 250, is shown in Fig. 1 by the
1000
black solid line.
900
B. Hybrid-Electric Vehicle Example
800
o
u = 10 xf The energy consumption of hybrid-electric vehicles can
700
be described well using a quasi-static discrete-time model.
state variable x [-]

600 The modeling follows the ideas described in [14], [15].

500 uo = 5 Essentially, the model contains the battery state-of-charge as
400
the only state variable. In a nutshell, the combustion engine is
modeled using an affine Willans approximation, the electric
300
x0 uo = 0 infeasible motor is modeled using an electric-power map (derived from
200 detailed simulations), and the battery is modeled as a voltage
100 source together with a resistance in series. The vehicle model
0 includes air drag, rolling friction, and inertial forces. The
0 20 40 60 80 100 120 140 160 180 200
gearbox is modeled using a constant efficiency of 95%.
time t [days]
The hybrid vehicle considered in this study has a 20%
Fig. 1. The optimal control signal map, determined using dynamic hybridization as defined in [16].
programming, for the discrete-time Lotka-Volterra system. The optimal state The model equations can be summarized and described as
xk+1 = f (xk , uk , vk , ak , ik ) + xk ,
trajectory for x0 = 250 when using the map is shown as the solid black
line. (21)
where xk is the battery state-of-charge, uk is the torque split
optimal control problem (15)–(20) is then solved with the factor, vk is the vehicle speed, ak is the vehicle acceleration,
dpm function using: and ik is the gear number. The model assumes isothermal
conditions of the components, no extra fuel consumption
% create grid during the startup of the combustion engine, and no energy
[Link]{1} = 201; losses during gear shifting. A constant auxiliary electric
[Link]{1}.lo = 0;
[Link]{1}.hi = 1000;
power demand of 350 W is used in the model.
[Link]{1} = 21; Since the drive cycle is assumed to be known in advance
[Link]{1}.lo = 0; the particular driving speed vk , acceleration ak and gear
[Link]{1}.hi = 10;
number ik at instance k can be included in the model
% set initial state function to form the time-variant model:
xk+1 = fk (xk , uk ) + xk ,
grd.X0{1} = 250;
k = 0, 1, . . . , N − 1. (22)
% set final state constraints
[Link]{1}.hi = 1000; The optimization problem of minimizing the total fuel mass
[Link]{1}.lo = 750; consumed
N −1
% define problem J = ∑ ∆mf (uk , k) ⋅ Ts (23)
[Link] = 1/5; k=0
prb.N = 200*1/[Link] + 1; for the hybrid vehicle over a given drive cycle, here the
% set options Japanese 10-15 driving cycle (J1015), can be stated as the
options = dpm(); discrete-time optimal control problem:
[Link] = 1;
N −1
min ∑ ∆mf (uk , k)
[Link] = 1;
[Link] = 1200; (24)
uk ∈Uk k=0
[Link] = 1;
s.t.
xk+1 = fk (xk , uk ) + xk
[res dyn] = dpm(@fishery,[],grd,prb,options);
(25)
x0 = 0.55 (26)
The output of the DP algorithm is an optimal control signal xN = 0.55 (27)
xk ∈ [0.4, 0.7]
map, specifying the optimal control signal at each time step k
and at each state xk ∈ Xk . The optimal control signal map for (28)
the Lotka-Volterra system is shown in Fig. 1. It shows that the 660
N= +1 (29)
optimal control is ”not fishing” u = 0 if the fish population Ts
is small x < 500, ”moderate fishing” u = 5 if the population where ∆mf ⋅ Ts is the fuel mass consumption at each time
is x = 500 and ”full fishing” u = 10 if the population is step. The time step in this example is Ts = 1 s. The optimal
large x > 500. Toward the end of the problem, one must control problem (24)–(29) is solved using DP. Figure 2
stop fishing as late as possible, such that the population shows the resulting optimal control map [Link]{1,:} and
reaches the specified minimum final size of xf,min = 750. state trajectory res.X{1} when using the dpm-function as
The resulting optimal state trajectory, i.e., the fish population described below.

1629
Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 19,2025 at 10:47:09 UTC from IEEE Xplore. Restrictions apply.
 two examples. This dpm function together with the model
functions introduced in this paper can be downloaded at [7].
0.65 The computational time3 required for backward calculation
0.6
for the two examples, without using the boundary line, is
state variable x [-]

shown in Fig. 3. It shows that the function evaluates 600000

5 0.55 points/s for the fishery problem and 200000 points/s for the
0.5 HEV problem. This is due to the more complex model in
the HEV problem. Future work includes the attempt of a
0.45
possible extension of the boundary line method to more
0.4 general problems, and the support of simple discrete-time
0.35
Simulink models. The main task for the near future will be
0 100 200 300 400 500 600 to optimize the memory requirements of the function.
time t [s]
VI. ACKNOWLEDGMENTS
Maanxd stliel

Torq
Rec
Stfeasib

Th e
In

Elec
Thanks to Daniel Ambühl and Moritz Oetiker for their
ha

rma

trica
u
rechl

e as
constructive inputs during the development of the code and
rge

l
arge

sist
for testing and evaluating the dpm function.
Fig. 2. Optimal input map obtained using the DP algorithm for a full R EFERENCES
parallel hybrid-electric vehicle driving the Japanese 10-15 drive cycle. The
[1] R. E. Bellman, Dynamic programming. Princeton - N.J.: Princeton
black curve shows the optimal state-of-charge trajectory when the battery
University Press, 1957.
is 55% charged at the start.
[2] D. Bertsekas, Dynamic programming and optimal control, 3rd ed.
Belmont, Massachusetts: Athena Scientific, 2005.
0.8
[points x106 /s]

[3] R. Luus, Iterative dynamic programming, ser. Monographs and surveys

HEV in pure and applied mathematics. Boca Raton: Chapman & Hall/CRC,
0.6 Fishery 2000, vol. 110.
0.4 [4] R. E. Bellman and E. S. Lee, “History and development of dynamic
programming,” IEEE Control Systems Magazine, vol. 4, no. 4, pp.
N ⋅N x⋅N u

0.2 24–28, 1984.

tcpu

[5] M. Back, S. Terwen, and V. Krebs, “Predictive powertrain control

0 for hybrid electrical vehicles,” in IFAC Symposium on Advances in
0 5 10 15 20 25 30
Automotive Control, Salerno, Italy, April 2004, pp. 451–457.
N ⋅ N x ⋅ N u [points x106 ]
[6] J. Pu and C. Yin, “Optimal control of fuel economy in parallel hybrid
electric vehicles,” Journal of Automobile Engineering, vol. 221, pp.
Fig. 3. The computational cost for the two examples. The values are given 1097–1106, 2007.
in calculated grid points per second as a function of the total number of [7] O. Sundström and L. Guzzella, “DPM-function,” Institute
grid points. for Dynamic Systems and Control, Department of
Mechanical and Process Engineering, ETH Zurich, 2009,
[Link]
% create grid [8] O. Sundström, D. Ambühl, and L. Guzzella, “On implementation of
[Link]{1}=61; [Link]{1}.hi=0.7; [Link]{1}.lo=0.4; dynamic programming for optimal control problems with final state
[Link]{1}=21; [Link]{1}.hi=1; [Link]{1}.lo=-1; constraints,” Oil & Gas Science and Technology - Revue de l’IFP,
% set initial state 2009, Accepted for publication.
grd.X0{1} = 0.55; [9] M. Schaefer, “Some aspects of the dynamics of populations important
% final state constraints to the management of the commercial marine fisheries,” Bulletin of
[Link]{1}.hi = 0.55; Mathematical Biology, vol. 53, pp. 253–279, 1991, Reprinted from the
[Link]{1}.lo = 0.55; Bulletin of the Inter-American Tropical Tuna Commission, 1(2):27–
% define problem 56, 1954.
prb.W{1} = speed_vector; % (661 elements) [10] H. Mosbech, “Optimal control of hybrid vehicle,” in International
prb.W{2} = acceleration_vector; % (661 elements) Symposium on Automotive Technology & Automation, vol. 2. Turin,
prb.W{3} = gearnumber_vector; % (661 elements) Italy: Automotive Automation Ltd, 1980, pp. 303–310.
[Link] = 1; [11] C.-C. Lin, H. Peng, J. W. Grizzle, and J.-M. Kang, “Power manage-
prb.N = 660*1/[Link] + 1; ment strategy for a parallel hybrid electric truck,” IEEE Transactions
% set options on Control Systems Technology, vol. 11, no. 6, pp. 839–849, 2003.
options = dpm(); [12] A. Sciarretta, M. Back, and L. Guzzella, “Optimal control of parallel
[Link] = 1; hybrid electric vehicles,” IEEE Transactions on Control Systems
[Link] = 1; Technology, vol. 12, no. 3, pp. 352–363, 2004.
[Link] = 1000; [13] A. Sciarretta and L. Guzzella, “Control of hybrid electric vehicles,”
[Link] = 5; IEEE Control Systems Magazine, vol. 27, no. 2, pp. 60–70, 2007.
[Link] = 'c'; [14] L. Guzzella and A. Sciarretta, Vehicle propulsion systems introduction
[Link] = 0; to modeling and optimization, 2nd ed. Berlin: Springer, 2007.
[15] L. Guzzella and C. H. Onder, Modelling and control of internal
[res dyn] = dpm(@hev,par,grd,prb,options); combustion engine systems. Berlin: Springer, 2004.
[16] O. Sundström, L. Guzzella, and P. Soltic, “Optimal hybridization in
two parallel hybrid electric vehicles using dynamic programming,”
V. CONCLUSIONS AND FUTURE WORKS in 17th IFAC World Congress, ser. Proc. of the 17th IFAC World
In this paper a Matlab function is introduced that effi- Congress, Seoul, Korea, 2008.
ciently solves deterministic DP problems. The syntax and
3
the main features of the function are highlighted using Calculated on a 32-bit Intel Pentium D 2.8GHz with 2.0 GB RAM.

1630
Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 19,2025 at 10:47:09 UTC from IEEE Xplore. Restrictions apply.