remarkRemark hypothesisHypothesis claimClaim Controlling Nutrition and
Propulsion Force in Running RacesC. Cook, S. Lenhart, W. Hager, and G. Chen

# Optimally Controlling Nutrition and Propulsion Force in a Long Distance
Running Race††thanks: Submitted to the editors 06/29/22. Support by the US
National Science Foundation under grant 2031213 and by the US Office of Naval
Research under grant N00014-22-1-2397 are gratefully acknowledged.

Cameron Cook RTI Health Solutions, Research Triangle Park, NC (,
https://www.rtihs.org/research-team/cameron-cook<EMAIL_ADDRESS>Suzanne
Lenhart Department of Mathematics, University of Tennessee, Knoxville, TN ().
<EMAIL_ADDRESS>William Hager Department of Mathematics, University of
Florida, Gainesville, FL<EMAIL_ADDRESS>Guoxun Chen Department of
Nutrition, University of Tennessee, Knoxville, TN<EMAIL_ADDRESS>
###### Abstract

Runners competing in races are looking to optimize their performance. In this
paper, a runner’s performance in a race, such as a marathon, is formulated as
an optimal control problem where the controls are: the nutrition intake
throughout the race and the propulsion force of the runner. As nutrition is an
integral part of successfully running long distance races, it needs to be
included in models of running strategies. We formulate a system of ordinary
differential equations to represent the velocity, fat energy, glycogen energy,
and nutrition for a runner competing in a long-distance race. The energy
compartments represent the energy sources available in the runner’s body. We
allocate the energy source from which the runner draws, based on how fast the
runner is moving. The food consumed during the race is a source term for the
nutrition differential equation. With our model, we are investigating
strategies to manage the nutrition and propulsion force in order to minimize
the running time in a fixed distance race. This requires the solution of a
nontrivial singular control problem. Our results confirm the belief that the
most effective way to run a race is to run approximately the same pace the
entire race without letting one’s energies hit zero.

###### keywords:

Optimization, Differential Equation Models, Running, Nutrition, Bioenergetics

34B60, 49N90, 92B05

##  1 Introduction

Running is one of the most popular forms of exercise. There are more than
$275,000$ road races per year in the United States [26]. At the marathon
distance alone, there are over $500,000$ people a year in the United States
who choose to race (www.runningtheusa.com). Some people run for fun while
others choose to seriously compete, attempting to run the shortest time for
the fixed distance race. Outside of sprints, pacing oneself to run the best
possible race is crucial. If the runner starts the race too slow, they may not
finish in as fast a time as expected. On the other-hand, if the runner goes
out too fast, they may find themselves running out of energy, struggling to
even finish the race. Knowing what pace strategy is best for each individual
runner is the hard part, as it depends on several values that are individual
to each runner, something this work aims to answer. Looking at data collected
from timing mats in the 2015 Boston Marathon, Figure 1 plots how the runners
placed versus how their pace differed from their first 10k and their overall
pace. Amongst this elite field, one can see that on average, the runners who
placed lower had much higher pace variation than runners who finished in a
shorter time (better placing).

Figure 1: The above data comes from the 2015 Boston Marathon. We graphed
finish placement on the x-axis and pace difference on the y-axis. Pace
difference is the difference between the runners 5km pace (min/mile) and the
runners overall pace (min/mile) for the race. The horizontal line represents a
difference of 0 between the two paces.

While running the fastest race possible has been a goal for runners, it has
also been a research topic for scientists for many years. Keller’s model
developed in the early 1970’s, was the first to cast running an optimal race
as an optimal control problem and has since been adapted and expanded by
others such as Aftalion and Bonnans [1], Woodside [31], and Pitcher [22]. All
of the results from these works show the importance of pacing, but lack
attention to the different energy systems. Also, no runner model to date,
considers in-race nutrition as an energy source available to the runner. In
order to not run out of energy (and run the best they can) during a long
distance run, a runner must consume food. There is a whole market of products
for runners to use in order to deliver necessary fuel to the body quickly.

Our goal of this work was to build a more dynamic runner model that better
represents the body’s energy systems, including energy allocation of fat and
carbohydrate energy dependent on velocity, as well as in-race nutrition. We
also wanted to minmize the time it takes to run different length races (in
particular the marathon distance race) by optimally choosing velocity and in-
race nutrition consumption profiles through control of propulsion force and
nutrition input. We first discuss the physics, biochemistry, and other factors
that can limit or enhance a runner’s performance before presenting our model,
and optimization techniques. Results and various simulations are presented,
followed by discussion.

##  2 Background for Building Our Model

###  2.1 Background and groundwork from past models

There are many different factors that coaches and runners have to think about
in preparing for and running a race. They have to consider all the training
needed to get the runner in the best possible condition for the race as well
as all the in-race components. It takes time to train not just the muscles,
but your body’s energy pathways and metabolism. All the cumulative knowledge
to develop the best training plan and best racing plan began with
understanding how the body transfers stored energy into mechanical work.
Originally, thermodynamics views exercise as the heat to mechanical energy
transfer. This idea was applied to interpretations of the two known energy
transfers: aerobic and anaerobic; however, biological energy transfer, or
bioenergetics, was discovered to be better suited to describe the anaerobic
and aerobic energy transfers due to chemical nature of the exchange.

Keller [15] considered Newton’s Second Law and oxygen supply in his model,
which was the first model of its kind, describing running races, and treats
energy as available oxygen in the muscles per unit mass. His first equation is
the equation of motion:

(1) $\frac{dV}{dt}=f-\frac{V}{\tau}$

where $t$ is time, $V(t)$ is the instantaneous velocity, $f(t)$ is the
propulsive force per unit mass , $\frac{V}{\tau}$ is a resistive force per
unit mass, and there is a constraint on the force, $0\leq f\leq f_{max}$. His
second equation governs energy, as the oxygen balance equation for $e_{an}$:

(2) $\frac{de_{an}}{dt}=\overline{\sigma}-fV$

where $\overline{\sigma}$ is the oxygen breathing and circulation rate in
excess of that supplied in non-running state with the constraint
$e_{an}(t)\geq 0$ for all $0\leq t\leq T$ [15]. Solving the optimization
problem and using only one energy type, Keller obtained results very close to
world records. For example in the 100m race, Keller’s theoretical results
garnered a time of 10.07 seconds while the world record at the time was 9.9
seconds, a percent error of $1.7\%$ [15]. As both anaerobic and aerobic
processes happen in the human body, modelers like Aftalion and Bonnans [1], as
well as others [31, 3]) included both processes in their model. Woodside [31]
added a fatigue factor for a long-distance runner model, but there is no
consideration of what measures runners can take to combat this. We aim to
address this and highlight the need to include 2 separate energies in a runner
model.

###  2.2 Quantifiable Values that drive Performance

When it comes to predicting one’s running race performance, there are many
“calculator” tools that have been developed because the best indicator for a
runner’s future race performance is their past race performances, as is
supported by decades of evidence in research [25]. Although there are new
predictors being developed using machine learning that need lots of data, most
runners have not run the number of races needed for such a robust data set.
Thus, we chose to focus on other physical indicators that we can model through
a system of equations that describe the physical and biochemical systems of a
runner.

Two quantifiable values that describe an athletes fitness are: VO2 max, and
VLamax. Both play a role in how an athlete expends energy; and more precisely,
what metabolic system they are primarily using at a given exertion level.
These values are comparisons of your aerobic and anaerobic systems. The
relative strengths of these systems determine how well an athlete will perform
at particular activities and can be used for good feedback in training.

Which metabolism is being used by the body is dependent on what percentage of
one’s current VO2 maximum (max) they are using, where VO2 max is the maximum
amount of oxygen you can utilize during exercise. The more fit you are, the
higher your VO2 max will be. Note that VO2 max is measured in milliliters of
oxygen consumed in one minute, per kilogram of body weight (mL/kg/min). This
value ususally ranges between 20-60 with professional athletes holding values
as high as 90 [28]. One can certainly improve their VO2 max through purposeful
training; although factors such as gender and genetics also play a role as the
size of one’s heart and their body fat percentage determine the body’s
capacity to carry oxygen through the blood. There are devices, formulas. and
VO2 max calculators [17, 29] that predict your VO2 max based on recent races
and current level of activity.

In preparing for long distance races, runners are attempting to improve their
VO2 max, as to allow the body to work at higher levels before needing to use
the anaerobic system as the energy pathway. This means running at faster
velocities without expending as much energy. When walking people are between
$15\%$ and $30\%$ of their VO2 max, solely using their aerobic process, with
the percentage increasing with exertion. When running a long-distance event,
the runner mainly uses aerobic respiration and aims to stay at their aerobic
threshold for the majority of the race. One’s aerobic threshold is at about
60% VO2 max and is a level of intensity that can be sustained for a long time,
where they are not using their anaerobic metabolism. Runners should only use
their anaerobic metabolism for a small portion of the race, as operating at
such an intensity cannot be sustained for very long. The anaerobic system
begins to be used between $75\%$ and $85\%$ of one’s VO2 max. During this
time, the aerobic system is still used, but at a lower rate as the percentage
of VO2 max increases up until $100\%$. At that point, the anaerobic system is
the main energy system with the aerobic system unable to be a significant
contributor as the individual is functioning at a level where they have used
up all available oxygen and need an energy system that does not require
oxygen.

A related, more tangible term to know is your velocity at your VO2 max, called
VVO2max. It is perhaps a more valuable piece of data to know than just knowing
your VO2max as it combines your aerobic capacity and running economy (or
running mechanics), into one comparable number [4]. To test for a runner’s
VVO2max, on a treadmill, their pace is increased, and the administrator
observes their VO2 increases linearly until it eventually starts to plateau.
When the VO2 plateaus, even as velocity is increased, no more oxygen is being
taken in and the runner has reached their VVO2max. So, VVO2max is the minimum
running velocity with peak oxygen uptake and it’s the speed at which there is
the lowest contribution from the anaerobic metabolism, while also having
maximum oxygen uptake.

Before reaching VO2 max, runners hit their lactate threshold. Lactate is
produced in anaerobic metabolism in the muscle whereas carbon dioxide and
water are created in aerobic metabolism. When the muscle lactate accumulates
before its removal, the intensity has to be lessened as the same level of
intensity is unsustainable [27]. Hitting one’s lactate threshold happens
between 70-85$\%$ of their VO2 max. This threshold is a crossover point of the
two energy processes where the body begins using anaerobic respiration for
energy. When the anaerobic metabolism becomes the main energy process, the
blood oxygen level is not enough to keep up with the energy demand.

A runner’s VLa or lactate capacity, is the body’s anaerobic power, or maximum
ability to produce lactate (G Hillson, personal communication, April 21,
2021). The higher the VLa, the worse the runner is at clearing lactate near
threshold. Marathon runners want to have a low VLa so that they can use more
fat for energy and spare their carbohydrates (G Hillson, personal
communication, April 21, 2021).

While VO2 max is a good indicator of fitness, a runner’s VLa max is what sets
the professional runner apart from other runners. When looking for athletes
for Nike’s sub 2 project [12], they only looked at athletes who had a VO2 max
greater than 65, but at that point VO2 max doesn’t tell the whole story. A
runner with a VO2 max of 65 or one with a VO2 max of 80, could have identical
optimal races, depending on their VLa max. It doesn’t matter if a runner has a
high VO2 if they aren’t able to access all of that VO2 (G Hillson, personal
communication, April 21, 2021). In conclusion, a long-distance runner wants to
train their body to have a high VO2max (or VVO2max) and a low VLamax.

###  2.3 Understanding the Body’s Energy Sources

The anaerobic metaboli-
sm provides an athlete with quick energy, through the creation of glucose, but
only 2 ATP (energy) molecules are obtained; wheras the aerobic metabolism
creates 38 ATP molecules through a more length process. The aerobic and
anaerobic systems occur in separate cellular compartments (mitochondria and
cytoplasm respectively) and often at different rates, involve different
reactants and products. Not only is the allocation of the two separate energy
processes of interest, but also which fuel is being utilized. Glucose and
fatty acids provide most of the fuel required for energy production in
skeletal muscles during aerobic exercise whereas glucose is the main source of
energy in anaerobic exercise. The body has significantly more energy available
in the form of fat, but the rate of using this energy form can not be
increased at high exercise intensities when the anaerobic metabolism is the
main mechanism. Thus, the body is mainly able to use fatty acids as an energy
source at low levels of intensity [30]. When you are not exercising,
approximately 30 percent of your energy comes from glycogen and 70 percent
from fat stores [30]. These percentages shift when intensity increases, as
does the number of calories being burned. Glucose is preferred as it is
readily available and quick metabolized, but is limited.

Figure 2 is a diagram showing fuel utilization between fat and glycogen as a
function of percent VO2 max. In Figure 2, at low percent of VO2 max, the fuel
utilization is low, and the percentage contribution from fat is substantial.
With the rise of VO2 max, when running faster, the rate of fuel utilization
increases and the percent contribution of fat decreases.

Figure 2: Fuel Utilization: allocation of calories per minute per kilogram
used as a function of percent VO2 max (adapted from [30]

###  2.4 Why We need Nutrition during the Race

When running races such as the marathon one must consider pacing strategies.
The runner must not only consider oxygen availability, but how much energy the
body has stored. To optimally run these races, runners are mainly using the
aerobic system, but their intensity with corresponding percentage of the VO2
max is high enough for the body to use glucose as energy, causing them to
still burn through the limited supply of glycogen. Runners of different
levels, masses, running at relatively different VO2 values, burn through this
glycogen at different rates; however, it is commonly accepted that on average
runners burn just over 100 kilocalories per mile, meaning their stores will be
depleted after about 20 miles or between 1.5 hours and hours [30, 20]. When
these stores are depleted, the term ‘bonking’ is commonly used, where a runner
must slow down significantly or forced to stop running altogether and walk. As
only one person to date has completed the marathon in under 2 hours, this is a
major problem for long distance runners that must combat. In the 1960’s, work
by [5] confirmed that blood glucose concentrations were linked to fatigue and
that eating hard candies during a race prevented weakness and fatigue during a
race [20]. Carbohydrates are digested in the small intestine and converted
into glucose. Glucose is stored mainly in the muscles and the liver as
glycogen, a chain of multiple glucose residues, but is also available for
immediate use if necessary [30]. The body can store about 600 grams of
glycogen, with 500 grams stored in the muscles and 100 grams stored in the
liver [30], totaling 2400 Kilocalories of glycogen stores in the body.

During a long distance running race when glycogen stores are depleted,
ingested or exogenous carbohydrates are quick sources of energy for the
muscles, available once absorbed by the muscles from the blood. Runners
typically consume gels, which are a concentrated dose of glucose throughout
long distances races. Most runners solely consume gels, that contain
approximately 25 grams of carbohydrates (100 Kilocalorie). Taking in nutrition
during the race allows a runner to move longer before glycogen stores are
depleted and they’ve reached some anaerobic energy threshold where they must
walk. It is known that the muscles absorb plasma glucose at a maximal rate of
1g/min-1.7 g/min [30, 14, 13] depending on the sugar mixture. This means that
while it may only take 3-5 minutes for some of the carbohydrates from a 100
calorie intake to reach your muscles, it can take approximately $25$ minutes
for all of the carbohydrates from the package to be absorbed. While there is
no limit to how many carbohydrates a runner can ingest, too many carbohydrates
consumed in a short time will result in digestive discomfort, forcing the
runner to slow down.

##  3 Our Model

We build on the current models by adding in some novel terms as well as
completely reformulating the energy equation. Our overall goal is to determine
the best nutrition and pacing strategy to use when running a marathon (or
other long distance race) in order to finish in the shortest amount of time.
Our first objective, is to develop a runner model that takes into account fuel
allocation depending on percentage VO2 max, fuel intake (in-race nutrition),
as well as the force which is applied to the ground by the runner, by using a
system of differential equations. Next, we would like to cast this as an
optimal control problem to determine the optimal velocity and in race
nutrition consumption profiles through control of propulsion force and
nutrition input.

Our system of ordinary differential equations for: $V,\text{Velocity,}\ \
E_{F},\text{Fat Energy,}\\\ E_{G},\text{Glycogen Energy,}\ \
N,\text{Nutrition}$ ,where $V$ will be measured in meters/min, $E_{F}$ and
$E_{G}$ in KJ/Kg, $N$ in KJ, and our control force $f(t)$ in meters/min2, are
given by:

(3) $\displaystyle\frac{dV}{dt}=\ $ $\displaystyle f(t)-\frac{V}{\tau}$ (4)
$\displaystyle\frac{dE_{G}}{dt}=\ $ $\displaystyle c_{3}j(N)-af(t)Vglyc(V)$
(5) $\displaystyle\frac{dE_{F}}{dt}=\ $ $\displaystyle-af(t)V(1-glyc(V))$ (6)
$\displaystyle\frac{dN}{dt}=\ $ $\displaystyle s(t)-dN-j(N)$

where $a=\frac{1}{1000}$. Figure 3 is a diagram of our system, including the
velocity compartment. Velocity is not directly connected to any of the other
compartments via energy transfer, but its impact on the two energy
compartments can be seen in equations (4) and (5). Our nutrition intake
strategy is the fuel source function, $s$, entering in the nutrition
compartment.

Figure 3: Energy-Nutrition-Velocity Diagram

In equation (3), there is a propulsive force per unit mass $f(t)$, and
$|\frac{V(t)}{\tau}|$ is a resistive force per unit mass. The initial
condition for this equation is: $V(0)=0$, as we start from rest, and we have a
constraint on the force, $0\leq f(t)\leq f_{max}$. This is the same equation
that Hill, Keller, and others used for velocity [9, 16, 31, 3, 1].

As energy is available from two different processes in the body, anaerobic and
aerobic metabolisms, we consider two energy sources. In aerobic metabolism,
the body uses energy from both fat and glycogen, while in anaerobic metabolism
only energy from glycogen is used, but at a faster rate. Thus, our total
energy, $E$, can be written as $E=E_{F}+E_{G}$. The human body has large
stores of fat for use as energy, but the body prefers to use glycogen when
moving at higher rates; thus, a higher percentage of fuel usage comes from fat
at lower velocities. As runners typically take gels that are mainly sugar (not
a significant source of fat), we assume there is no input into the fat energy
compartment. The body expends energy at a rate, $f(t)V$ (work done). As we are
allocating the energy usage between the two sources, our fat energy
differential equation (5), and glycogen energy differential equation (4) have
glycolitic function, $glyc(V)$, a fuel allocation function of one’s glycogen
energy through the anaerobic and aerobic pathways dependent on the quotient of
one’s instantaneous velocity, $V$, and their velocity, $VVO2max$, at 100$\%$
VO2max, written $\frac{V}{VVO2max}$. Thus, $1-glyc(V)$ accounts for the fuel
allocation function of one’s fat energy through the aerobic pathway dependent
on velocity compared to one’s velocity at 100$\%$ VO2max.

It is difficult to know one’s current percentage VO2 max; however, one could
know a priori their velocity at 100 percent VO2 max. Billat and Koralsztein
[4] showed that one’s VO2 max and their velocity at VO2 max are linearly
related. As velocity is more easily computable we will use this relationship
in our model in determining the allocation of energy at different velocities
compared with one’s velocity at VO2 max. To approximate $glyc(V)$, we first
described it it as a piecewise continuous function on intervals of
$\frac{V}{vvmax}$ based on 2 from [30] and then approximated that function by
a function that smooths the points where the derivatives do no exist
(described in detail later). Equations (4) and (5) have a convex combination
of this function and the fat allocation function of $1-glyc(V)$, as the total
work rate must equal $f(t)V$. We obtained the piecewise graph for $glyc(V)$
from Figure 2 as well as from the literature [30]. Figure 2 gives a good
estimate of how the body uses fat energy versus glycogen energy. We assume
that the anaerobic pathway is not used until after $60\%$ VO2 max, and that
the energy from fat linearly decreases to 0 by $100\%$ VO2 max at which point
the anaerobic system, and thus the glycogen compartment, is solely used.
Figure 4 is a graph of three different possible $glyc(V)$ functions that we
obtained from our Figure 2.

Recall, that while what percentage of one’s VO2 max at which they are running
is explicitly in the $glyc(V)$ function, VO2 max is not the only value that
effects our fuel utilization. Lactate capacity (VLa), also impacts one’s fuel
utilization as described earlier. We account for differences in VLa, by
considering different structures for our fuel utilization function, $glyc(V)$.
The function is still dependent on percent VO2max, but varies in fuel
utilization at particular percent VO2 max values. Figure 4 shows three
different $glyc(V)$ functions corresponding to 3 different lactate capacities.
We use these 3 different glyc structures to represent runners who have a low
VLa (good), an average VLa (avg), and a high VLa (bad). The runners with a low
VLa (good) are able to run at a faster pace than those with a higher VLa
without accumulating as much lactate in their muscles. The VLa is a feature
that sets the best professionals apart from one another.

In equation (4) for glycogen energy there is a source term with $j(N)$ that
comes from the nutrition differential equation (6). Our $j(N)$ term is a
nutrition consumption function increasing energy available in the muscles in
the form of glycogen, but at a bounded rate: $j(N)=c_{4}N$ with rate constant
$c_{4}=\frac{1}{m}$, the inverse of the runners mass. Our initial conditions
for our two energy equations are $E_{G}=144$ Kilojoules per unit mass,
assuming the runner has full glycogen stores, and $E_{F}(0)=3439$ Kilojoules
per unit mass, dependent on the runner’s body fat percentage. Both of the
energies must stay non-negative, and dictate the choices of $f(t)$ and the
corresponding V. Thus, in particular, if $E_{G}=0$, $f(t)=0$. This is not
optimal and would be avoided in an optimal race until the very end.

Figure 4: Three different glyc functions corresponding to runners with a good,
average, or bad VLa.

Further, in our nutrition energy differential equation, (6), $s(t)$ is a
source term from nutrition input such as a gel. Each gel is roughly 100
Kilocalories (4.18 Kilojoules). The addition of the source term slows the rate
at which $E_{G}$ decerases. The $-dN$ term in the nutrition differential
equation represents the nutrition used for basic bodily function, not
available to the muscles for energy usage. The initial condition for this
equation is $N(0)=0$ as we assume there is nothing in this compartment when
the race starts. Model parameters can be found in the results section in Table
2.

##  4 Optimal Control

Our second goal is to determine the best strategy to use when running a long
distance race in order to finish in the shortest amount of time. Casting this
as an optimal control problem, we could think of this as a problem of
maximizing the distance over a fixed time interval or equivalently (proved in
[1]), as a problem of minimizing time over a fixed distance. We have two
controls in our problem: propulsive force, $f(t)$, and fuel intake $s(t)$. We
will determine the optimal control $f(t)$ for each of the intake strategies
and then optimize over those intake strategies. If we chose to solve the
minimum time problem. an extra isoperimetric constraint must be considered.
Thus, we choose to solve the maximum distance problem.

###  4.1 Optimal Control Problem Setup

For the maximum distance problem, assuming we have a nutrition strategy,
$s_{i}$, we maximize the objective functional:

(7) $J_{i}(f)=\int_{0}^{T}V(t)dt$

for $s_{i}\in S$ (a finite set of nutrition strategies), where

$S=\\{s_{1},s_{2},...,s_{N}\in L^{2}(0,T)|0\leq s_{i}(t)\leq
s_{max},a.e.,i=1,2,...,N\\}$

we determine an optimal control depending on $s_{i}$, and then optimize over
our set $S$. We obtain continuous velocity and energy profiles corresponding
to our control force and then optimize over discrete $s_{i}$. We solve the
following with fixed $T$:

(8) $\max_{S}\max_{A_{i}}J(s_{i},f)=\max_{S}\max_{A_{i}}\int_{0}^{T}V(t)dt$

with bounded controls: $0\leq f\leq f_{max}$, $0\leq s_{i}(t)\leq s(t)_{max}$,
and control set

$U_{i}=s_{i}\times\\{f\in L^{2}(0,T)\,|\,0\leq f(t)\leq f_{max},a.e.\\}$

subject to state constraints: $0\leq E_{F}(t)$ , $0\leq E_{G}(t)$, for
$t\in[0,T]$
initial conditions: $V(0)=0$, $E_{F}(0)=K$, $E_{G}(0)=144$, $N(0)=0$,
and state equations (3)-(6). Our admissible control set for $s_{i}\in S$ is

$A_{i}=\\{(s_{i},f)\in U_{i}\,|\,E_{F}(t)\geq 0,E_{G}(t)\geq 0,\
\text{for}\,t\geq 0\\}$

Note that $V(t)$ and $N(t)$ are non-negative, from (3,6).

###  4.2 Existence of the optimal control

######  Theorem 4.1.

Given a nutrition strategy $s_{i}$, there exist an optimal control $f^{*}$
that solves the optimal control problem of maximizing (7) with corresponding
state equations (3)-(6) over admissible controls,

$\mathcal{U}=\\{f\in L^{2}(0,T)\,|\,0\leq f(t)\leq f_{max},a.e.,\,E_{F}(t)\geq
0,\ E_{G}(t)\geq 0,\ \text{for}\ t\geq 0\\}.$

######  Proof 4.2.

First, the admissible control set is nonempty since the control $f=0$
satisfies the conditions. Our objective functional, $J(f)$ is uniformly
bounded due to bounds on the states and controls. There exists a maximizing
sequence of control functions denoted by ${f^{n}}$ such that

$\lim_{n\rightarrow\infty}J(f^{n})=\sup_{f\in\mathcal{U}}J(f),$

and a corresponding sequence of states $V^{n},E_{F}^{n},E_{G}^{n},N^{n}$. Note
that the control sequence and the state sequences are uniformly bounded. By
the Banach-Alaoglu theorem [8], there exists a control $f^{*}$ such that on a
subsequence $f^{n}\rightharpoonup f^{*}$ weakly in $L^{2}(0,T)$. From the
state differential equations we see that the derivative sequences are also
uniformly bounded. Thus, $V^{n},E_{F}^{n},E_{G}^{n},N^{n}$ are uniformly
Lipschitz and thus equicontinuous. Therefore by the Arzela-Ascoli theorem, on
a subsequence, we have uniform convergence of our corresponding state
functions

$V^{n}\rightarrow V^{*},E_{F}^{n}\rightarrow E_{F}^{*},E_{G}^{n}\rightarrow
E_{G}^{*},N^{n}\rightarrow N^{*}\ \ .$

We need to show that $f^{*}$ is an optimal control and
$V^{*},E_{F}^{*},E_{G}^{*},N^{*}$ are corresponding optimal states. To show
those states correspond to control to $f^{*}$, we illustrate this by showing
the $E$ differential equation:

$E_{F}^{n}(t)-E_{F}(0)=\frac{-1}{1000}\int_{0}^{t}f^{n}V^{n}(1-glyc(V^{n}))ds$

and

$\displaystyle|\int_{0}^{t}\frac{f^{n}V^{n}}{1000}(1-glyc(V^{n}))ds-\int_{0}^{t}\frac{f^{*}V^{*}}{1000}(1-glyc(V^{*}))ds|$
$\displaystyle\leq\ \ $
$\displaystyle||f^{n}||_{2}||V^{n}-V^{*}||_{2}+|\int_{0}^{t}V^{*}(f^{n}-f^{*})ds|\
\ \ \ (\text{by Holder's Inequality})$

We have that all 4 states correspond to the control $f^{*}$ . Next, $f^{*}$ is
an optimal control since,

$\displaystyle\lim_{n\rightarrow\infty}J(f^{n})$
$\displaystyle=\lim_{n\rightarrow\infty}\int_{0}^{T}V^{n}dt=\int_{0}^{T}V^{*}dt=J(f^{*})\
.$

Now that we’ve proven an optimal control exists, we begin solving the optimal
control problem using Pontryagin’s Maximum Principle (PMP) [23] by determining
the set of necessary conditions, using the Hamiltonian with state constraints.
We solve the maximum distance problem, obtaining an $f^{*}$ for each nutrition
strategy, $s_{i}$, and maximize over our set of finite strategies $S$.

###  4.3 Necessary Conditions for the Maximum Distance Problem

We determine the necessary conditions by using PMP [23] to compute the
Lagrangian, which is the Hamiltonian with state constraints:

$\displaystyle\mathcal{L}$
$\displaystyle=V+\lambda_{1}(f(t)-\frac{V}{\tau})+\lambda_{2}(-af(t)V(1-glyc(V)))+\lambda_{3}(c_{3}j(N)-af(t)Vglyc(V))$
$\displaystyle+\lambda_{4}(s(t)-dN-j(N))+\eta E_{G}\ ,$

where the penalty function $\eta(t)\geq 0$ is a Lagrangian multiplier appended
to $E_{G}$ as a state constraint penalty. Note that we do not have a state
constraint penalty for $E_{F}$ as $E_{F}$ will not get close to $0$ over the
time interval. The function $\eta$ satisfies: $\eta\equiv 0$ where
$E_{G}^{*}>0$, and $\eta\geq 0$ otherwise when the state constraint is tight
(i.e. $E_{G}=0)$ Since we seek to maximize $\mathcal{L}$ with respect to
$f(t)$, a state variable violating the constraint would decrease $\mathcal{L}$
and not be optimal as $\eta E_{G}<0$.

Also, PMP[23] gives the existence of adjoint functions , $\lambda_{1}$,
$\lambda_{2}$,$\lambda_{3}$ and $\lambda_{4}$, satisfying system (9) (11)
appending the state equations to the Hamiltonian.

(9) $\displaystyle\lambda_{1}^{\prime}=-\frac{\partial\mathcal{L}}{\partial
V}$
$\displaystyle=-(1-\frac{\lambda_{1}}{\tau}-\lambda_{2}af(t)+\lambda_{2}af(t)V\frac{\partial
glyc}{\partial V}+\lambda_{2}f(t)aglyc(V)-$
$\displaystyle\lambda_{3}af(t)V\frac{\partial glyc}{\partial
V}-\lambda_{3}af(t)glyc(V))$ (10)
$\displaystyle\lambda_{2}^{\prime}=-\frac{\partial\mathcal{L}}{\partial
E_{F}}$ $\displaystyle=0\ \ ,\ \
\lambda_{3}^{\prime}=-\frac{\partial\mathcal{L}}{\partial E_{G}}=-\eta$ (11)
$\displaystyle\lambda_{4}^{\prime}=-\frac{\partial\mathcal{L}}{\partial N}$
$\displaystyle=-(c_{3}c_{4}\lambda_{3}-\lambda_{4}d-\lambda_{4}c_{4})$

with transversality conditions $\lambda_{1}(T)=0$, $\lambda_{2}(T)=0$,
$\lambda_{3}(T)=0$, $\lambda_{4}(T)=0$

For $t\in[0,T]$, due to it’s initial condition, $E_{F}$, stays positive for
any choices of $f(t)$ in the admissible control set and thus $\lambda_{2}=0$
over $[0,T]$. Since our objective functional and state equations are linear in
the control, we consider the different signs of the derivative of
$\mathcal{L}$ with respect to $f(t)$, which is the switching function
$\phi(t)$.

(12) $\phi(t)=\frac{\partial\mathcal{L}}{\partial
f}=\lambda_{1}-\lambda_{3}aVglyc(V)$

At time $t$, for $\phi(t)<0$ we will have a maximum when $f^{*}(t)=0$, whereas
for $\phi(t)>0$ we will have a maximum when $f^{*}(t)=f_{max}$. However, PMP
[23] does not tell us what happens when $\phi(t)=0$. More information with
respect to control $f$ than $\frac{\partial\mathcal{L}}{\partial f}=0$ is
needed when $\phi(t)=0$ on a subinterval. If this happens only for for finite
number of time points, we would have a bang-bang control and those points
would represent the switching times. On the other hand, if $\phi(t)=0$ on a
subinterval of time, then we would have a singular control.

For optimal control problems with state constraints there is the possibility
of boundary singular sub-arcs as well as interior singular sub-arcs. Boundary
singular sub-arcs occur when the constraint is tight, meaning when
$E_{G}(t)=0$. If $E_{G}(t)=0$, that implies $E_{G}^{\prime}=0$, and solving
(4) to obtain $f_{s,b}$, provided that $V(glyc(V))\neq 0$. Note that as
$glyc(V)$ accounts for the percent of glycogen energy being used and there is
always a percentage of fuel usage coming from glycogen, this function will
never be $0$. Due to the structure of the differential equation for $V$, as
soon as $f$ is positive, $V$ is also positive. Hence, $V(glyc(V))\neq 0$ when
the boundary singular arc is active.

If the constraint is not tight, then we have the possibility of an interior
singular sub-arc that we can obtain by differentiating the switching function
with respect to time twice, noting that on an interior singular sub-arc
$\eta(t)=0$, and solving for $f_{singular}$.
When the constraint $E_{G}(t)$ is not tight, $\eta(t)=0$; however, when it is
tight we solve for it using the fact that over the small time interval where
the glycogen energy constraint is tight, the switching function, $\phi(t)$, is
$0$ on a boundary subarc and thus also $\phi^{\prime}(t)=0$. We are then able
to solve this equation for $\eta(t)$ to get our complete characterization of
$\eta(t)$.

For singular subarcs, one incurs an additional necessary condition known as
the generalized Legendre-Clebsch (GLC) condition [24]. The necessary GLC
condition for the existence of a singular subarc for an optimal control
problem, in which the control appears in $\phi^{\prime\prime}(t)$ is

$\frac{\partial}{\partial f}\frac{d^{2}}{dt^{2}}\frac{\partial}{\partial
f}(\mathcal{L})\ =\ \frac{\partial}{\partial f}\frac{d^{2}\phi}{dt^{2}}\ \ \
\geq\ \ 0,$

a condition one can verify for this problem. For our last necessary condition
we have that the Transversality Conditions:
$\lambda_{1}(T)=\lambda_{2}(T)=\lambda_{4}(T)=0\ .$ From the above, our
optimal control, force, has the following structure:

$\displaystyle f(t)=\left\\{\begin{array}[]{ccc}0&for&\phi(t)<0\\\ &&\\\
f_{singular}(t)&for&\phi(t)=0\,\ E_{G}\neq 0\\\ &&\\\
f_{s,b}(t)&for&\phi(t)=0\,\ E_{G}=0\\\ &&\\\ f_{max}&for&\phi(t)>0\\\
\end{array}\right.$

This $f(t)$ represent the force profile of which a runner should adhere to in
order to run the optimal race. We have shown that our problem meets all of the
necessary conditions for existence of a singular interior sub-arc, including
the GLC condition, which suggests that the optimal trajectory could include a
singular component. To summarize the control trajectory, we know that the
optimal control begins with a maximum force sub-arc, that for a race with
large enough $T$, an optimal control must be comprised of more than just a
maximum force arc, that a singular boundary force sub- arc exists (and likely
at the end of the time interval), that a singular interior sub-arc is likely,
and that intuitivly, it is unlikely for a zero force arc to exist, as that
would not be optimal. We believe that for our system, the control is singular
for the majority of the event, and is comprised of a maximum force sub-arc,
followed by a singular interior sub-arc, and finishing with a singular
boundary sub-arc, as has been shown to be the trajectory in simpler systems in
[15, 22, 1, 31]. Unfortunately one is unable to easily obtain an exact
solution structure with switching times. As the problem is linear in the
control, and believed to be largely singular, along with having state
costraints and the nutrition pulse function, it is challenging to solve even
using adjoints and penalty multipliers. Thus, we chose to discretize the
system in order to implement obtaining an approximate solution to the optimal
control problem. In the next section we discuss our nutrition strategies. We
then outline our approximation method where we discretize the problem and add
a penalization technique to reduce the noise in the solution.

###  4.4 Nutrition Strategies

Not only are we concerned with what force $f$ value will maximize the distance
over a fixed time T, but also with which in race nutrition strategy, $s_{i}$
is optimal. As there are many philosophies within the running community about
how often you should take nutrition during a race, we test 15 different
nutrition strategies, $S=\\{s_{0},s_{1},s_{2},\dots,s_{15}\\},$ shown in Table
1.

During Eulid Kipchoge’s race to the beat the 2 hour marathon barrier, his
nutrition consisted of hydrogels by Maurteen. He consumed 100g Carbs per hour
(which would be 4 baked sweet potatoes in solid food for reference) [12],
altough there is no comment on his exact feeding regime. This amount of
carbohydrates is pretty high for a runner to tolerate during a marathon due to
the increased gastrointestinal movement in running. This number of
carbohydrates is more in line with what cyclist would consume during their
longer races. It is expected that at his level, his nutrition provided him a
1-3 percent boost in performance compared to him solely consuming water during
the race. This is in contrast to the practice most runners follow of taking
about 50 carbs per hour, an area of research that is continuously being
studied by nutrition companies. We have many different nutrition strategies
over which we are optimizing, two of which include simulating Kipchoge’s known
race intake, and a standard intake strategy used by average runners. We will
look at the percentage improvement from using the determined strategy versus
consuming no carbohydrates. As a reminder, we do not have stomach sensitivity
to food throughout the race in our model yet, which could skew the results in
favor of taking as much nutrition in as possible.

Table 1: Scenarios of nutrition intake during a marathon race $s_{0}$,$s_{1}$,$s_{2}$,$s_{3}$, $s_{4}$,$s_{5}$,$s_{6}$, $s_{7}$ |  100 Calorie gels spread evenly 0,1,2,3,4,5, 11, and 24 times throughout the race.
---|---
$s_{8}$ |  One 100 Calorie gel taken towards the beginning of the race.
$s_{9}$ |  One 100 Calorie gel taken towards the end of the race.
$s_{10}$ |  Four 100 Calorie gels: 2 taken early in the race, 2 taken toward the end of the race.
$s_{11}$ |  Four 200 Calorie gels spread evenly throughout the race.
$s_{12}$, $s_{13}$ |  250 Calorie gels spread evenly 2,4 times throughout the race
$s_{14}$ |  Ten 50 Calorie gels taken evenly throughout the race.
$s_{15}$ |  One 250 Calorie gel taken, 1 100 Calorie gel taken, 1 250 Calorie gel taken.

##  5 Methods

Over time, optimal control problems have been solved and approximated using
many different numerical techniques. There are many numerical examples in
mathematical biology that use the Forward Backwards sweep method [10, 18].
Programs such as GPOPS and PASA have been developed to handle particular types
of optimal control problems [11, 21]. MATLAB has a minimization tool, fmincon,
that is built to handle a variety of optimization problems. We tried to use
the forward backwards sweep method, the packages GPOPS and PASA, as well as
fmincon in the continuous setting; however, none of these methods were robust
enough to handle this particular problem. Thus, we discretized our problem and
used fmincon, inputting our differential equation system in through equality
constraints. We began the discretization of our optimal control problem by
partitioning our time interval, $[0,T]$ using M+1 equally spaces nodes,
$0=t_{0}<t_{1}<\cdots<t_{M}=T$. We used a left rectangular approximation for
the objective functional, obtaining the maximization problem:

(14) $\max_{f}J(f)=\max_{f}[\sum_{k=0}^{M-1}hV_{k}]$

where $h=\frac{T}{M}$ , with bounded controls: $0\leq f_{k}\leq f_{max}=36000$
meters/min2, and $f=(f_{0},\dots,f_{M-1})$, with $s=(s_{1},\dots,s_{15})\in
S$, and $0\leq s_{i}(t_{k})\leq s_{max}$ for all $0<k<M-1$ , and with initial
conditions: $V_{0}=0$, $E_{F,0}=K$, $E_{G,0}=144$, $N_{0}=0$.

Next we use a forward Euler approximation for the state equations and obtain:

(15) $\displaystyle V_{k+1}$ $\displaystyle=V_{k}+h(f_{k}-\frac{V_{k}}{\tau})$
(16) $\displaystyle E_{F,k+1}$
$\displaystyle=E_{F,k}+h(-af_{k}V_{k}(1-glyc(V_{k})))$ (17) $\displaystyle
E_{G,k+1}$ $\displaystyle=E_{G,k}+h(c_{3}j(N_{k})-af_{k}V_{k}(glyc(V_{k})))$
(18) $\displaystyle N_{k+1}$ $\displaystyle=N_{k}+h(s(t_{k})-dN_{k}-j(N_{k}))$

We also now have the discretized version of $J(N_{k})$: $j(N_{k})=c_{4}N_{k}\
.$ We discretized our glyc function, $glyc(V_{k})$, by using a spline
interpolator in MATLAB, described in the next section, and dropped in each
$i^{th}$ nutrition scenario, $s_{k}$, corresponding to the time points each
nutrition strategy designated. To optimize our system, while satisfying the
constraints, we use fmincon as the minimization solver on the discretized
system, using a left-rectangular integral approximation for the objective
function.

“Fmincon is a gradient-based method that is designed to work on problems where
the objective and constraint functions are both continuous and have continuous
first derivatives” [19], but due to our constraints and nutrition input, it
was helpful to write it as a discrete system. Fmincon is set to accept: the
objective function, a starting guess for the control vector, X0, inequality
linear constraints, equality linear constraints, lower and upper bounds for
the state variables, nonlinear inequality and equality constraints, as well as
options that allow the user to change the MATLAB settings for various features
or provide the system with more information. Our discretized objective
function was originally

$J=-\delta\sum_{i=1}^{T-1}V_{k}$

where $V_{k}$ occurs in our $X$ vector as entries $X(3n+1:4n)$ and
$\delta=\frac{T}{M}$ where $T$ is the length of the race and $M$ is the number
of mesh points, including $t=0$, but we will modify this objective function to
include minimizing the total variation in $f$. In order to modify the
objective function to include minimizing the total variation, we added 2
additional state variables, $\iota$ and $\zeta$ [2], such that the variation
in 2 time points in $f$ is written as

$f(i+1)-f(i)=\zeta(i)-\iota(i)\ .$

Our starting guess vector, $X0$, that gives an initial placeholder for our
state vector was: $X0=[\overbrace{36000\ \dots 36000}^{T}\ \overbrace{0\ \dots
0}^{6T-2}].$ Note that each state vector variable is a column vector of length
$(T\times 1)$ and the two vectors used for variation penalization, $\zeta$ and
$\iota$ are of length $(T-1)\times 1$ due to their structure. Our initial
condition vector was: $[f(1)\ E_{F}(1)\ E_{G}(1)\ V(1)\ N(1)\ \zeta(1)\
\iota(1)]=[36000\ 3439\ 144\ 0\ 0\ 0\ 0]\ .$ As we do not have any inequality
constraints and our equality linear constraints were written as:
$A_{eq}x=b_{eq}$, where $Aeq$ is a $n\times n$ matrix, $x$ is our solution
column vector of length $n\times 1$, and $beq$, of length $n\times 1$ is the
righthand side of our equality linear constraints.

As the velocity and nutrition equations are linear, we input those as our
linear equality constraints. Our 2 energy equations are non-linear and thus
are input as non-linear equality constraints. One important feature for our
problem is our energy constraints, $E_{G}(i)\geq 0$ and $E_{F}(i)\geq 0$ for
every $i\in[0,T]$. With our problem set up in its current format, we are able
to simply input our lower bounds for our fat and glycogen energies as 0. We
input the runners physical force, velocity, and energy capacities as
appropriate upper bounds. We also give our variables appropriate initial
conditions for their starting out fat and glycogen energies, 0 velocity, and 0
in race nutrition.

In our glycogen energy equation we have a function we call “glyc” and its
counter part “1-glyc” in our fat energy equation. We approximated these from a
fuel utilization Figure 4 adapted from one in Stipanuk and Caudill’s textbook
[30], by creating 2 vectors of points where one vector represents the percent
VO2 max and the other represent the percent glycogen used. To smooth out our
glyc function, we used a function in MATLAB called spline. Spline is a cubic
piecewise polynomial interpolator that is continuous and twice differentiable
everywhere. It takes the 2 vector of points you supply it with and creates
$n-1$ cubic polynomials that connect at the supplied $x$ vector given.

One common issue in control problems where the control has a singular sub-arc
is variation. Ding and Lenhart [7] as well as Caponigro et al. [6] both
penalized their original objective functional to regularize the chattering. In
order to obtain a reasonable trajectory for our runner, we chose to penalize
our objective functional by adding a term to our objective functional that
bounded the total variation [2]. Adding this penalty, not only amends the
objective function to include penalizing for variation, it adds a 5th
constraint to our problem. By adding this penalty we are decreasing the total
variation in the solution trajectory. We can express the total variation,
$\mathcal{V}(f)$, in the control [2] as:

$\mathcal{V}(f)=\sup\sum_{i=0}^{N-1}|f(i+1)-f(i)|$

where the supremum is being taken over all possible partitions of our time
interval. Due to the difficulty of differentiating absolute value functions we
decompose the absolute value term, using 2 new $T-1$ vectors $\zeta(i)$ and
$\iota(i)$, whose entries are non-negative. Each entry of $\zeta(i)$ and
$\iota(i)$ will be defined as

$f(i+1)-f(i)=\zeta(i)-\iota(i)$

We can think of this decomposition as satisfying the following 2 conditions
[2]:

If $f(i+1)-f(i)>0$, then $\zeta(i)=f(i+1)-f(i)$ and $\iota(i)=0$;
If $f(i+1)-f(i)\leq 0$, then $\zeta(i)=0$ and $\iota(i)=-(f(i+1)-f(i))$,
where $\zeta(i)\geq 0$, $\iota(i)\geq 0$ and either $\zeta(i)=0$ or
$\iota(i)=0$

Appending this penalty to the objective function we have:

$J=-\delta\sum_{i=1}^{T-1}V(i)+p\sum_{i=1}^{T-2}(\zeta(i)+\iota(i))\ ,$

and resulted in the extra linear equality constraint

$ceq5(i)=-f(i+1)+f(i)+\zeta(i)-\iota(i)\ .$

The coefficient $p$ of our bounded variation scales the degree to which we
penalize variation. If $p$ is large, there is a more emphasis being placed on
minimizing the variation of the control variable, force.

##  6 Results and Discussion

###  6.1 World record optimal results

We found that using a direct discrete optimization without adjoint functions,
fmincon, was the most efficient approximation method and captured all the
dynamics. Table 2 shows the constants and coefficients with their scientific
meaning and their units. Five parameters that can be chosen depending on the
individual are: $E_{G}(0)$ (Initial glycogen energy), mass, VLa type,
nutrition uptake rate,$c_{4}$,and VVO2 max, shown at the end of the table. We
first chose to simulate the optimal race of the current world record holder.
In the attempt to break the 2 hour marathon barrier, experts collaborated to
optimize a set of runners VLa max, VO2 max, general running economy, nutrition
plan, as well as pick the perfect race course for these runners. The parameter
values for our simulation can be found in Table 2. Note that the 3 different
VLa types correspond to runners with 3 different VLamax’s. These particular
parameter values were picked to reflect appropriate values for a professional
(world record breaking) runner. Plots of the state and control variables are
seen in Figure 5.

In Figure 5 we see that the runner’s velocity rapidly increases from 0 to
reach an approximately constant velocity of $357$ m/min, that they can
maintain for the entirety of the race. The plot of the runner’s propulsion
force, similar to the runner’s velocity plot, shows the runner’s force quickly
increasing from 0 m/min2 to a constant force of $2.14\cdot 10^{4}$ m/min2. The
runner’s fat energy, $E_{F}$, decreases linearly, from $3439$ KJ/Kg to $3390$,
nowhere near the energy constraint, $E_{F}\geq 0$. The runner’s glycogen
energy decreases from its initial value of $E_{G}(0)=144$ KJ/Kg to its final
value of approximately $0$ KJ/Kg. We see bumps in this subplot corresponding
to the nutrition fueling. In this simulation we assumed that the runner took
four 200 calorie gels at $t=20$, $46$, $71$, $97$, which can be seen in the
nutrition subplot in Figure 5, and the runner took 120 minutes to complete the
distance.

Table 2: Parameters Used for World Record Holder Parameter | Value | Unit | Meaning
---|---|---|---
$T$ | $120$ | minutes | length of race
$M$ | $T+1$ | minutes | number mesh points including $t=0$
$\tau$ | $1/60$ | min | internal resistant force constant
$d$ | $0.005$ | min-1 | loss of nutrition
$sm$ | $1/3600$ | (seconds)-2 | seconds to minutes conversion
$c_{3}$ | $1/m$ | kg-1 | mass conversion constant
$p$ | $0.5$ | unitless | variation penalty coefficient
$a$ | $1/1000$ | kilojoules/joules | unit conversion constant
$\delta$ | $120/121$ | T/M | discretization parameter
$m$ | $55$ | kilogram | runner mass
$E_{G}(0)$ | $144$ | kilojoules/kg | initial energy at start of race
$c_{4}$ | 1/6 | min-1 | nutrition uptake rate
VVO2max | 402 | meters/min | Velocity at 100% VO2max
VLa type | good | unitless | good, average, bad in glyc function
Figure 5: World record marathon simulation, with four 200 calorie gels

During the world record breaking race, the runner did not consume gels in the
same manner that we implement, but did consume an equivalent amount of energy.
The VVO2max was set to be $402$ meters per minute. It is expected that world
class runners have a VVO2 max in this range. This VVO2max is approximately
equivalent to a runner running $4$min/mile. We see that the distance achieved
in our results is slightly longer than that of a marathon. The fastest
marathon ever run thus far is 1:59:40. In the simulation from our model, the
runner runs 120 minutes in 42.5 km, which is equivalent to the runner running
a 1:59:09 marathon. This is only a $0.4$ percent difference between the
current world record and our results. This difference in average pace would be
the runner running 4:33 minutes per mile instead of 4:34 minutes per mile.
Note, that although our results are extremely close to the current world
record, parameters could be changed to address a runner having an even better
VLa, or better running economy. We also didn’t explicitly have “shoe type” in
our model but a better shoe could result in better running economy which would
result in a better VVO2 max. It is believed that the world record time will
continue to get faster over time.

###  6.2 Results: Varying runner dependent parameters

A runner’s in-race nutrition, VLamax, and their VVO2max are the driving
factors behind a runners performance, and thus it is important to see how
changing these impacts the results. We first describe the results we obtained
from simulations of our system with each of the 15 different nutrition
strategies. The parameters used for these simulations are all the same ones
used in Table 2, except, we changed $T=135$, which also changes
$\delta=\frac{135}{136}$.

Recall our 15 different nutrition strategies, as well as taking no nutrition
labeled $s_{0}$, from Table 1. Table 3 shows the distance completed with each
of these strategies. There are several conclusions we can draw from these
results, which we will present after looking at graphs of several of the
nutrition strategies. First we compare a strategies of taken 0 nutrition
versus five 100 cal gels. From Table 3 we can see that the runner is able to
run between $0.5$ and $0.7$ kilometers further per each additional gel, which
translates to approximately $1.5\%$ increase in distance per gel. This means
that by taking five 100 calorie gels instead of taking no nutrition, the
runner has a $7.75\%$ improvement in performance. In Figure 6 and Figure 7 we
see similar strategies in the force and velocity subplots; however the
proplsion force and therefore the velocity that the runner who took 5 gels is
able to maintain is significantly higher than that of the runner who took no
nutrition.

Figure 6: Marathon simulation with 3 states and optimal force with no gels  Figure 7: Marathon simulation with 4 states and optimal force with five 100 calorie gels Table 3: Total Distance Achieved with Different Nutrition Strategies  Strategy | Distance (km) | Strategy | Distance (km)
---|---|---|---
$s_{0}$ | 40.0 | $s_{8}$ | 40.7
$s_{1}$ | 40.7 | $s_{9}$ | 40.5
$s_{2}$ | 41.4 | $s_{10}$ | 42.5
$s_{3}$ | 42.0 | $s_{11}$ | 45.1
$s_{4}$ | 42.6 | $s_{12}$ | 43.2
$s_{5}$ | 43.1 | $s_{13}$ | 45.6
$s_{6}$ | 46.0 | $s_{14}$ | 43.1
$s_{7}$ | 52.9 | $s_{15}$ | 43.7

We also checked if there would be any differences in performance if the runner
took in the same amount of calories, but distributed differently. In $s_{5}$,
where the runner took in 500 calories by taking five 100 calorie gels the
distance achieved was 43.1 km. In $s_{12}$ where the runner takes in two 250
calorie gels the distance achieved was 43.2 km. Lastly, in $s_{14}$ where the
runner consumes ten 50 calorie gels, the runner traveled 43.1 km. So,
regardless of the way in which the runner consumes the same amount of
calories, they travel relatively the same distance. Perhaps this would be
different if there were a mechanism in the model that discouraged eating too
much at once or eating later in the race (we attempt to remedy this in an
extension of our model). Next, we compare strategy $s_{5}$, where the runner
consumes a total of 500 calories, and strategy $s_{8}$, where the runner
consumes a total of 800 calories.

A runner taking five 100 calorie gels throughout the race is around the
average number of gels used by your typical marathon runner. Some of the
professionals are taking up to 100 carbohydrates an hour, which translates to
400 calories an hour. A runner finishing the marathon in about 2 hours would
take in 800 calories, which is simulated with strategy $s_{8}$, in Figure 5.
In the simulation the runner is able to run $1.5$km further, an improvement of
$3.5\%$. In both simulations we see similar trajectories; however, the runner
in simulation $s_{8}$ is able to run at a higher average speed throughout the
race, yielding the larger distance ran in the same amount of time. These
results were expected as increasing the energy the body is able to use, should
result in greater distance traveled. This shows the importance in runners
being able to take in as many calories during a race as possible. It is
possible to increase the amount of calories one can take in during a race. One
can do this by taking nutrition that is easier to digest, taking water with
the nutrition so that the salinity balance of the gut is kept, and practicing
taking in nutrition during long training runs.

Finally we take a look at the strategies, $s_{6}$, $s_{7}$ where the runner
takes eleven and twenty-four 100 calorie gels throughout the race
respectively. It is important to mention that these scenarios are unlikely to
occur as a runner would struggle to consume this many gels without getting
sick. We completed simulations with nutrition strategies $s_{6}$ and $s_{7}$
to see if the model would organically capture the effect of a runner taking in
too much nutrition. As can be seen in the results from Table 3, this effect
was not captured and the simulations simply show an increase in energy as more
gels are added, and no negative side effects. In the simulation where the
runners consumes 11 gels, they are able to run 46.0 km, while in 24 gel
simulation the runner travels 52.9 kms. Simply put, with a fixed time race, we
see an increase in distance traveled for every gel that is taken. Taking in
eleven and certainly twenty-four 100 calorie gels throughout a race would most
certainly upset the runners stomach and all of those carbohydrates would
struggle to make it to the muscles due to imbalances in the stomach. The body
does not handle the amount of accumulated energy in the nutrition compartment
without a negative reaction. In future work we will handle this issue.

Next, we chose to vary the VLa types. We use the same parameters from Table 2,
and vary the glyc function.

Table 4: Total Distance Achieved by Runners with Different VLa Types across Two Nutrient Strategies VLa type | Nutrition | Distance (km)
---|---|---
good | 0 calories taken | 41.7
average | 0 calories taken | 40.0
bad | 0 calories taken | 37.8
good | Four $100$ calorie gels | 44.1
average | Four $100$ calorie gels | 42.6
bad | Four $100$ calorie gels | 40.5

Table 4 shows the distances achieved using different VLa types across two
different nutrition strategies. The VLa types labeled as “good”, “average”,
and “bad” correspond to the glyc functions shown in Figure 4. For these runs,
we tested the distance for the case where the runner takes in no nutrition
during the race, as well as the runner taking in four $100$ calorie gels. In
the three different simulations with varying VLa’s, the runner consumes four
100 calorie gels during the race. In these results, the importance of runners
training their body’s to have lower maximum rate production of lactate for a
long distance runner is evident. The runner with the good VLa type ( i.e. the
runner whose $glyc$ function was made to use fewer carbohydrates than fats
than the runners with average or bad VLa’s) is able to run much further in the
same amount of time. There is an $8.9\%$-$10.3\%$ improvement in the runner
with a good VLa as compared with a bad VLa depending on if the runners took in
nutrition as well. We see the same race structure regardless of the VLa type,
with both simulations showing the runners ending with approximately $0$
glycogen energy, but similarly to the simulation where the runner takes more
nutrition, the runner with the higher VLa is able to maintain a higher force.
The runner with a good VLa can maintain a force of $1.97\times 10^{4}$ m/min2
throughout the race, while the runner with a bad VLa was only able to maintain
a level force of approximately $1.81\times 10^{4}$ m/min2. Some of these
results can be seen in Table 5

Lastly, we consider varying the runner’s VV02max. As this term applies to all
levels of runners, we vary many of the parameters including: the length of
race, $T$, the runners mass, $m$, the initial glycogen energy, $E_{G}(0)$, the
runners velocity at $100\%$ VO2 max, VVO2max, and their VLa type. We use many
of the same parameters from Table 2, but change the weight of the runner, $m$,
length of race, $T$, as well as the VVO2max of the runners. Note that $c_{3}$
and $M$ also vary, but do so due to their relationship with $m$ and $T$
respectively. Runners of different levels will finish the marathon in
different lengths of time, thus the need to vary $T$. Also, runners of
different abilities can have significantly different VVO2max values and VLa
types, and therefore change depending on the individual. Table 5 shows the
distances achieved for our different levels of runners with 4 lengths of time.
We completed these simulations for marathon runners finishing in different
lengths of time using fmincon optimization.

For runners who are running for 155 minutes and 180 minutes with VVO2max’s of
320 and 250 respectively, we see realistic results. The runners run at around
$85\%$ of their VO2 max that is suggested by experts. We continue to see the
impact of taking in-race nutrition compared without nutrition. The runner who
takes 180 minutes to finish the race and has an average VLa has a $4\%$
improvement when they take nutrition versus when they do not. One issue with
these results is with the “slowest runner”, since they are running at a pretty
high percentage of their VO2 max. It is unlikely that a runner would be able
to do that for this length of race. This is more likely due to the fact that
our model does not penalize a runner for running too fast, as long as they
have enough glycogen energy. Recall, that our only constraints are that the
energies have to be above $0$, and that there are some upper bounds on forces
and velocities. We will amend our current model to work for shorter distance
races to address this issue of runners only being constrained by total amount
of glycogen energy.

Table 5: Total Distance Achieved with Different Levels of Runners and Four Lengths of Time  Race Time (min) |  VVO2 max (m/min) |  weight (Kg) |  VLa type |  Nutrition (Cals) |  $E_{G}(0)$ (KJ/Kg) |  Dist. (Km)
---|---|---|---|---|---|---
$155$ |  320 |  73 |  avg |  4 100cal gels |  150 |  42.9
$155$ |  320 |  73 |  avg |  0 100cal gels |  150 |  40.1
$155$ |  320 |  73 |  good |  4 100cal gels |  150 |  44.4
$155$ |  320 |  73 |  bad |  4 100cal gels |  150 |  41.6
$180$ |  250 |  80 |  avg |  0 100cal gels |  140 |  41.0
$180$ |  250 |  80 |  avg |  4 100cal gels |  140 |  42.7
$215$ |  200 |  80 |  bad |  0 100cal gels |  140 |  41.6
$210$ |  200 |  80 |  avg |  4 100cal gels |  144 |  44.3

##  7 Conclusions and Future Work

To understand the impact of nutrition we formulated our model such that in-
race nutrition could be accounted for in the energy differential equations. We
created a marathon model described by velocity, fat energy, glycogen energy,
and in-race nutrition differential equations. This marathon model better
represented the body’s energy systems and included nutrition pulses throughout
the race. After formulating an optimal control problem with this system of
differential equations, and with the goal to maximize the distance achieved
over a fixed time interval, we proved the existence of an optimal control. We
determined a discretized approximate solution to the associated optimal
control problem using MATLAB’s fmincon optimization software, that included a
bounded variation penalty on our control force. We then optimized over a
finite set of in-race nutrition strategies, obtaining a solution to our model
for each strategy.

We analyzed the complicated role that in-race nutrition input has on a
marathon runner as well as energy allocation dependent on the percentage of
one’s VVO2max during running. Expressing the dynamics of the body’s energy
throughout a running race as 2 differential equations, representing fat energy
and glycogen energy, with usage determined by an allocation function dependent
on the ratio of one’s velocity compared to their VVO2max, is a valuable
contribution to the field of runner models. Prior to our model, all energy
dynamics during running races had been described based on available oxygen per
unit mass does not allow for one to include in-race nutrition energy. As a
runner has a much smaller storage supply of glycogen than fat, and the body
prefers to use the glycogen at high levels of exertion, the in-race nutrition
consists of essentially only carbohydrates. These factors required us to split
up the energy into the two compartments and to use the allocation function.
While we obtained an optimal solution for an expert runner, we also
individualized our model such that it can be used to determine a race strategy
for any type of runner as long as they know their weight, VVO2max and VLamax
type.

To obtain our approximate solution to the optimal control problem of
maximizing distance for a fixed time, we atte discretized our problem, and
then optimized, using Matlab’s optimizing tool, fmincon. We included a
penalization for variation in the objective function. The variation
penalization was extremely important and the solution structure was very
sensitive to this penalization. With too small of a variation penalty
coefficient, the runners purposive force was often erratic and non-optimal. We
found that the optimal trajectory for marathon runners was to run at a steady
state velocity for the majority of the race and take in at least four 100
calorie supplements spread evenly throughout the event. Over a variety of
runner’s individual parameters, this racing strategy was optimal. A runner can
improve their performance with consistent and targeted training to improve
their VV02max and their VLamax, (which correlates to having a better $glyc$
function) which can in turn drastically improve the their performance
potential. When we approximated a solution to our problem for a runner with
the parameters of the current world record holder in the marathon, our
approximation was quite good with an error of $0.4\%$.

Our contributions to the field are: formulating a system that includes energy
allocation dependent on the ratio of their velocity and their velocity at
VO2max, adding in-race nutrition to a runner model, and approximating a
solution to the optimal control problem, amended to include a penalty bounding
the variation of the control. As we have found that taking more carbohydrates
throughout the race improves performance, it is natural to believe that one
should take in as many carbohydrates as possible; however, this is unrealistic
as the stomach’s digestion slows through a running race and also struggles to
tolerate high volume of carbohydrates in the form of glucose and fructose. Our
initial marathon model does not account for this fact and our results showed
that a runner’s performance increased linearly with each gel they consumed. To
account for this limitation requires additional work, which will be presented
in our next paper.

Upon executing different scenarios for our marathon model, we realized that
our model was not appropriate for races between a sprint and a half marathon.
In these scenarios we obtained results that were much faster than has been
accomplished by humans thus far. These results came from only having energy
constraints in model, when in reality there is a difference in energy in the
body and energy that is available for use in the body, dependent on lactate
production. Our extended model with the effects of lactate will be presented
in our followup paper. In the future we would like to incorporate hydration
into our model. Hydration levels can be a huge factor for a runner. Low
hydration can negatively affect the rate at which the body uptakes the
nutrition. One possibility could be changing our nutrition uptake function to
represent low, medium, or high levels of hydration. Building this feature in
would be the next step in advising runners on optimal intake during running
races.

##  8 Bibliography

## References

  * [1] Aftalion, A. and Bonnans, J. F. (2014). Optimization of running strategies based on anaerobic energy and variations of velocity. SIAM Journal on Applied Mathematics, 74(5):1615-1636.
  * [2] Atkins, S., Aghaee, M., Martcheva, M., and Hager, W. (2020). Solving singular control problems in mathematical biology, using PASA, arXiv:2010.06744, 2020.
  * [3] Behncke, H. (1993). A mathematical model for the force and energetics in competitive running. Journal of Mathematical Biology, 31(8):853-878.
  * [4] Billat, L. V. and Koralsztein, J. P. (1996). Significance of the velocity at vo2max and time to exhaustion at this velocity. Sports Medicine, 22(2):90-108.
  * [5] Burke, L. M., Hawley, J. A., Wong, S. H. S., and Jeukendrup, A. E. (2011). Carbohydrates for training and competition. Journal of Sports Sciences, 29(sup1):S17-S27. PMID: 21660838. 30
  * [6] Caponigro, M., Ghezzi, R., Piccoli, B., and Trelat, E. (2018). Regularization of chattering phenomena via bounded variation controls. IEEE Transactions on Automatic Control, 63(7):2046-2060.
  * [7] Ding, W. and Lenhart, S. (2009). Optimal harvesting of a spatially explicit fishery model: Optimal harvesting. Natural Resource Modeling, 22(2):173-211.
  * [8] Friedman, A. (1982). Foundations of Modern Analysis. Dover Books on Mathematics Series. Dover.
  * [9] Furusawa, K., Hill, A. V., and Parkinson, J. L. (1927). The dynamics of “sprint” running. Proceedings of the Royal Society of London. Series B, Containing Papers of a Biological Character, 102(713):29-42.
  * [10] Hackbusch, W. (1978). A numerical method for solving parabolic equations with opposite orientations. Computing, 20(3):229-40.
  * [11] Hager, W. W. and Zhang, H. (2016). An active set algorithm for nonlinear optimization with polyhedral constraints. Science China. Mathematics, 59(8):1525-1542.
  * [12] Hutchinson, A. (2020). We now have the lab data on Nike’s breaking 2 runners. https://www.outsideonline.com/health/training-performance/nike-breaking2-runners-lab-data/
  * [13] Jentjens, R. L. P. G. and Jeukendrup, A. E. (2005). High rates of exogenous carbohydrate oxidation from a mixture of glucose and fructose ingested during prolonged cycling exercise. British Journal of Nutrition, 93(4):485-492. 205
  * [14] Jeukendrup, A. E. and Wallis, G. A. (2005). Measurement of substrate oxidation during exercise by means of gas exchange measurements. International Journal of Sports Medicine, 26(S 1):S28-S37.
  * [15] Keller, J. B. (1973). A theory of competitive running. Physics Today, 26(9):43-47.
  * [16] Keller, J. B. (1974). Optimal velocity in a race. The American Mathematical Monthly, 81(5):474-.
  * [17] Koutlianos, N., Dimitros, E., Metaxas, T., Cansiz, M., Deligiannis, A., and Kouidi, E. (2013). Indirect estimation of vo2max in athletes by acsm’s equation: valid or not? Hippokratia, 17(2):136-140.
  * [18] Lenhart, S. and Workman, J.T. (2007) Optimal Control Applied to Biological Models. Mathematical and Computational Biology Series, Chapman and Hall/CRC, London.
  * [19] MATLAB Optimization Toolbox (2020). Matlab optimization toolbox. The MathWorks, Natick, MA, USA.
  * [20] Murray, B. and Rosenbloom, C. (2018). Fundamentals of glycogen metabolism for coaches and athletes. Nutrition Reviews, 76(4):243-259.
  * [21] Patterson, M. A. and Rao, A. V. (2014). Gpops-ii: A MATLAB software for solving multiple-phase optimal control problems using Hp-adaptive gaussian quadrature collocation methods and sparse nonlinear programming. ACM Transactions on Mathematical Software, 41(1):1-37.
  * [22] Pitcher, A. B. (2009). Optimal strategies for a two-runner model of middle-distance running. SIAM Journal on Applied Mathematics, 70(4):1032-1046.
  * [23] Pontryagin , L. S., Boltyanskiii, V. G., Gamkrelize, R. V., and Mishchenko, E. F. (1962) The Mathematical Theory of Optimal Processes, Wiley.
  * [24] Robbins, H. M. (1967). A generalized Legendre-Clebsch condition for the singular cases of optimal control. IBM Journal of Research and Development, 11(4):361-372.
  * [25] Roberts, J. W. (2004). Predicting marathon performance using body mass variables and training data: A field prediction model. ProQuest Dissertations Publishing.
  * [26] www.runningtheusa.com (2021).
  * [27] Scott, C. (2005). Misconceptions about aerobic and anaerobic energy expenditure. Journal of the International Society of Sports Nutrition, 2(2):32-37.
  * [28] Shete, A. N., Bute, S. S., and Deshmukh, P. R. (2014). A study of vo2 max and body fat percentage in female athletes. Journal of Clinical and Diagnostic Research, 8(12):BC01-BC03.
  * [29] Stipanuk, M. H. and Caudill, M. A. (2005). Daniels’ Running Formula.
  * [30] Stipanuk, M. H. and Caudill, M. A. (2013). Biochemical, physiological, and molecular aspects of human nutrition. Elsevier/Saunders, St. Louis, Mo, 3rd ed. edition.
  * [31] Woodside, W. (1991). The optimal strategy for running a race (a mathematical model for world records from 50 m to 275 km). Mathematical and Computer Modeling, 15(10):1- 12.

## Acknowledgements

We would like to acknowledge the help from Greg Hillson, Bobby Holcombe, and
Summer Atkins.