O PTIMIZATION of human performance by accurately modeling fatigue has challenged athletes, coaches, and scientists. The rise in popularity of wearable sensors in physical activity tracking presents opportunities for modeling and optimizing performance as they alleviate the need for expensive laboratory equipment. Understanding fatigue dynamics can potentially help athletes train in a more efficient way and perform at their peak. Moreover, it can provide useful information to further enhance the performance of an athlete during a physical exercise. Fatigue due to prolonged exercise is defined as a decline in muscle performance which accompanies a sensation of tiredness [2] , [3] . Therefore, during physical exercise, fatigue prevents athletes from producing the required power. Several studies such as [4] [5] [6] have investigated fatigue in cycling and developed models for it based on the expenditure of anaerobic energy. Only a few studies investigate the recovery dynamics of this anaerobic energy [7] . Figure 1 illustrates a cyclist pedaling in a time trail on a mountain course. Even seasoned cyclists may find it challenging to pace themselves in such a course. By over exerting themselves too early or too late, they may not achieve their maximum potential. We hypothesize that cyclists can finish a time trial faster if they plan in advance in consideration of their Anaerobic Work Capacity (AW C) and upcoming road elevation. The pacing strategy can be formalized as an optimal 1 : An illustration of optimal pacing of a cyclist using upcoming elevation data, a dynamic model of cyclist fatigue and recovery, and a model of bicycle longitudinal dynamics. The illustration also shows a dynamic programming grid on two states of velocity and anaerobic energy for planning the optimal velocity or power. The optimal power is communicated and displayed to the cyclist in real-time. Drawn in https://www.icograms.com. control problem which requires dynamic models of fatigue and recovery and maximal power capacity of a cyclist as a function of their fatigue levels. While former members of our team had attempted optimal pacing based on a dynamic model of muscle fatigue and recovery in [8] and [9] , such lumped muscle models were complex and hard to experimentally verify and calibrate. Alternatively, one can consider fatigue as running out of AW C which could be captured with a single dynamic state and recovery is achieved when pedaling below a Critical Power (CP ). Such a model can be more easily verified and calibrated in laboratory experiments as our recent work in [1] and [10] shows. In fact, a model based on anaerobic energy expenditure has been used in [11] to optimize a 5 km cycling time-trial. However, a model for recovery was not considered in [11] perhaps because the road was assumed to be flat, not necessitating rest.
In this paper, the dynamic models for a cyclist's fatigue and recovery during anaerobic exercise presented in [1] are further validated and calibrated for six human subjects using 4 hours of experimental data per subject that we collected in the lab. We believe these models based on anaerobic work capacity are more practical to use in optimal pacing than those based on estimating muscular fatigue presented in [8] and [9] . This is mainly because the state of AW C can be estimated rather reliably open-loop and by integrating the cyclist's power output over time and without the need for invasive (e.g. blood lactate) measurements. Furthermore, we show in this paper that affine dependence of system's dynamics and constraints on control input (pedaling power) allows for valuable analytical insights into the nature of optimal pacing strategy when using Pontryagin's Minimum Principle. We resort to numerical solution via Dynamic Programming (DP) to resolve the switching between various modes of optimal strategy. We simulate one of the subjects on four time-trial courses of increasing level of difficulty including stage 13 of 2019 Tour De France.
Section II presents a literature review on muscle fatigue and recovery. Section III introduces our proposed mathematical models for anaerobic energy expenditure and recovery besides maximum power generation ability of cyclists. In Section IV the experimental protocol and results for seven cyclists are presented. After developing the models we formulate the optimal pacing problem in Section V. The problem is investigated analytically in Section VI using optimal control theory, and numerical solutions using dynamic programming are shown in Section VII, followed by conclusions.
Several studies have investigated fatigue, which has led to its many definitions. For example, in [?] authors define fatigue as a reduction in maximal capacity to generate force or power during exercise. Whereas, fatigue is defined as inability to produce the desired force or power resulting in impaired performance in [13] and [14] . The point of occurrence of fatigue is defined as the moment at which a drop from the desired power level is observed, which is also time-toexhaustion, thus making both terms indistinguishable [15] . To address this, in [16] and [17] fatigue is defined as a continuous process altering the neuromuscular functional state resulting in exhaustion and exercise termination. In general, there are two sources of fatigue: central fatigue and peripheral fatigue [18] . Central fatigue is defined as the failure of the Central Nervous System (CNS) in sending necessary commands to muscles to operate [19] . This factor is essentially important in high-intensity exercise [20] . Peripheral fatigue is caused at the muscle level. It can be induced because of the neuromuscular failure in muscles to comprehend and perform the commands coming from CNS, and deficiency of vital substances [20] .
The above mentioned studies illustrate the difficulty involved in arriving at a global definition of fatigue. The mechanism of fatigue leading to exhaustion are different for cycling, running, swimming, etc [21] . As we previously presented in [1] , we define muscle fatigue and recovery in cycling as follows:
• Fatigue: Expending energy from anaerobic metabolic systems by pedaling above a critical power which results in a decrease of maximum power generation ability. • Recovery: Recovering energy from anaerobic metabolic systems by pedaling below the critical power which results in an increase of maximum power generation ability. 
Measuring the effect of central fatigue on muscle performance is a challenge since the the matter is highly subjective [22] . Most research efforts to objectively measure central fatigue are focused on measurement of Maximum Voluntary Contraction (MVC) [23] [24] [25] . In voluntary contraction of a muscle, the generated force is proportional to the muscle electrical activation [26] . A standard way of measuring muscle activity is via Electromyography (EMG) tests. During the test the amount of electric potential produced in muscles can be measured. Although studies such as [27] and [28] have shown the accuracy of EMG tests in measuring maximal and submaximal voluntary contractions, there is evidence of underestimating the muscle activation at high force levels [18] . Therefore, to point to the goal of the current study, EMG cannot provide accurate data for modeling a cyclist's fatigue and recovery.
In addition to central fatigue there are several sites for peripheral fatigue. To get a better understanding of fatigue at muscle level, we should focus on muscle metabolic system. Like all of the natural processes in the world, muscle contraction needs a source of energy to be executed. The fuel that provides this energy is adenosine triphosphate (ATP). When one phosphate radical detaches from ATP, more than 7300 calories of energy are released to supply the energy needed for muscle contraction [29] . After this detachement, ATP converts to adenosine diphosphate (ADP). When a human muscle is fully rested, the amount of available ATP is sufficient to sustain maximal muscle power for only about 3 seconds, even in trained athletes [29] . Therefore, for any physical activity that lasts more than a few seconds, it is essential that new ATP be formed continuously.
There are three metabolic systems which provide the needed ATP: Aerobic system , Glycogen-lactic acid system and Phosphagen system. Table I compares the three systems, in terms of moles of ATP generation per minute and endurance time at maximal rates of power generation. Utilization of these systems during physical activity is based on the intensity of the activity.
Cycling can fall in all categories above depending on the intensity of the exercise. Many people cycle for fun and get around cities. In this case they may only use their aerobic system which provides them with low amount of power which they can hold for a very long time. However, during high intensity cycling such as in a time trial, the human body will use the other two sources besides the aerobic system to provide enough energy for muscle contraction. This is important when designing mathematical models that describe muscle power generation in cycling. Later on, we discuss a method to define a power limit below which the cyclists use their aerobic metabolic system that allows them to hold their power for long periods of time. However, there is limited energy to pedal above this power limit.
During aerobic exercise, muscles utilize the aerobic metabolic system to produce ATP. Oxygen plays a vital role in formation of ATP molecules. There are two major methods to measure oxygen during a fatiguing exercise. The first one is measuring the volume of oxygen intake in breathing (V O2 ) [30] . During a physical exercise,V O2 increases to provide the necessary oxygen needed to produce ATP in the muscle as suggested in [31] , [32] , and modeled in [33] . The experimental procedure to measureV O2 requires a number of laboratory equipment that cannot be used by a cyclist during everyday outdoor ride. The second method is directly measuring the amount of oxygenated hemoglobins at the local muscle (muscle oxygenation). When muscles are fully fresh, the percentage of oxygenated hemoglobins among the total number of hemoglobins is at its highest. During a fatiguing exercise, this percentage drops. Several studies show that Near Infrared Spectroscopy (NIRS) is a robust method to measure muscle oxygenation [34] [35] [36] . A few companies currently make wearable devices that enable the cyclists to monitor their muscle oxygenation in real-time [37] , [38] . We have shown in [1] a real-time measurement of muscle oxygenation during a set of experiments.
Another fatigue indicator is the amount of lactate produced in the muscle. Maximum Lactate at Steady State (MLSS) is the maximum maintainable blood lactate concentration without additional buildup through aerobic exercise [39] , [40] . During an anaerobic exercise, the rate of lactate production is higher than its dissipation [41] . As authors in [42] and [43] suggest, the amount of lactate accumulation can provide useful information about the fatigue level of muscles. The high lactate levels at muscle represent a lower ability of the muscle to generate force and power [13] , [44] , [45] . Traditionally, a common way to measure lactate has been taking blood samples or biopsies during an experiment [46] , [47] . These invasive methods are not in any way suitable for a real-time measurement and estimation of fatigue. Recently, authors in [48] developed a non-invasive method to measure blood lactate using electromagnetic wave sensors. Commercialization of such non-invasive in-situ measurement techniques will enable researchers to develop mathematical models that represent the relationship between blood lactate level and muscle power/force generation capacity.
In cycling, it is easier to measure power without elaborate laboratory equipment owing to the development of commercial grade power meters. However, determining exercise intensity using power is not straightforward as a threshold power is needed to classify exercise intensity. As stated in [49] and [50] , the Critical Power (CP ) can potentially be used to determine exercise intensity. It has been shown that CP is close to the power at which MLSS occurs according to [51] [52] [53] . This means, while pedaling at a power level above CP , a cyclist would expend energy from anaerobic energy sources. On the other hand, pedaling below CP helps the cyclist recover the Fig. 2 : The 3-minute-all-out test protocol. The average power at the last 30 seconds of the tests is considered to be CP , and the area between power plot and CP is equivalent to AW C. amount of energy burned using anaerobic metabolic system [51] [52] [53] .
The critical power concept was introduced by authors in [54] . They defined CP as the maximum power output that can be maintained indefinitely. In [55] authors showed that there is a limited amount of anaerobic energy for a cyclist to pedal above CP . This "tank" of energy is called Anaerobic Work Capacity (AW C). They suggest by pedaling at a certain power level above CP , a cyclist can hold that power for a limited amount of time before he or she runs out of anaerobic energy and cannot pedal above CP . This is known as the two parameter model. The relationship between critical power and anaerobic work capacity is often expressed as:
where P is instantaneous power in Watts, and t lim is time-toexhaustion in seconds. The experimental protocol designed to calculate CP and AW C in Equation (1) requires multiple lab visits for the test subjects according to [5] , [6] , [56] [57] [58] .
To avoid these multiple lab visits, authors in [59] developed a 3-minute-all-out test (3MT) as in Figure 2 . In this test, subjects sprint "all-out" for the entire 3 minutes. The value of CP is given by the average power output of the last 30 seconds. The value of AW C is the area between the power curve and CP . This test has been validated in [60] and [61] .
As discussed in Section II-C, AW C is a finite energy store for pedaling above CP . When a cyclist expends his or her AW C entirely, the maximum power that can be produced is CP . Let w be the amount of AWC remaining. The rate of change of w while expending energy above CP is given by the difference between the rider's power and CP . Recovery rate can be calculated similarly as assumed in [62] . However, it has been shown in [63] and [64] that recovery occurs at a slower rate than expenditure. For instance, if the cyclist is pedaling at a power level below the critical power, the amount of recovered energy will be less than the area between CP and the power curve. Thus, we propose to compute an adjusted recovery power using the amount of energy recovered given by,
where P adj is the adjusted recovering power, and T rec is the duration of recovery. Equation (2) can be rewritten to form an energy expenditure and recovery model as,
where P is the power output of the cyclist. Equation (3) represents a switching dynamic model of fatigue and recovery. To calculate the adjusted power we need a model that can provide P adj at any power level in the recovery interval (i.e. below CP ). One of the assumptions of this study is that the recovery rate of w does not depend on the duration of the recovery interval T rec . This assumption stems from the need for a causal energy recovery model. If it is assumed that w's rate of recovery depends on T rec , the amount of recovered energy during recovery interval will depend on the recovery duration in the future which does not seem plausible. In our previous study in [1] , we elaborated more on the reason behind making this assumption by using muscle oxygenation data.
To arrive at a relationship between P adj and the actual applied power by the rider, P adj calculated from the tests at each recovery power is plotted against the actual pedaling powers, and the following linear model provides a good fit to the data we present later in the paper,
where a and b are the model's constants. Another parameter which is affected by the expenditure of AW C is the cyclist's ability to generate an instantaneous maximum power at any point during cycling. This is illustrated in the 3MT where it is required for the subjects to pedal at their maximum power throughout the test. The power continuously decreases as a subject expends energy from AW C. Authors in [11] assumed a constant maximum power generation ability of 800 Watts regardless of the amount of energy remaining. If this assumption was to be valid, the 3MT profile would look completely different.
Using 3MT, a model of the maximum power vs. remaining energy w can be determined. The remaining energy at any instant during 3MT can be calculated by computing the remaining area between the power plot and the CP line. The remaining energy w can then be plotted against power values in 3MT. Maximum power of the cyclist is applied 2 to 3 seconds into the test as it can be seen in Figure 2 . Data points before the maximum power account for overcoming flywheel and aerobic inertia and are removed when fitting the model. The data is then fitted with a linear expression of the form,
where w is the remaining energy of the cyclist's work capacity, and α is the model's constant.
The experimental protocol comprised of three tests namely, (i) a ramp test to determineV O2 and Gas Exchange Threshold (GET), (ii) a 3-minute all-out test (3MT) to determine CP and AW C, and (iii) an interval cycling test to determine the recovery of AW C. The ramp test involves incrementally increasing the power until the subject is exhausted. From the ramp test, oxygen uptakeV O2 , defined as the volume of oxygen inhaled per minute per kilogram of body weight [30] , is determined. Furthermore, during the test, there is an abrupt change in the ratio of volume of CO 2 exhaled to the volume of O 2 inhaled, which is the point at which blood lactate concentration starts to increase. The power atV O2 (maximum oxygen uptake) and GET are recorded to aid in designing the 3MT test. The CP and AW C from the 3MT are also used to design the interval test for modeling the AW C recovery. The experimental protocol was approved by the Institutional Review Board of Clemson and Furman Universities. Besides power, muscle oxygenation and heart rate were also recorded during the tests. The interval test was developed using the definitions of fatigue and recovery form Section II-A to derive mathematical models for recovery of AW C as well as the maximum power generation capacity of cyclists at any instant in a ride. The protocol was developed under the following assumptions:
• The 3MT test provides reliable estimates of CP and AW C. • Exercise below CP utilizes the aerobic energy source, and results in recovery of AW C. • The recovery of AW C below CP happens at a slower rate than its expenditure above CP . • The rate of recovery depends on the recovery power and not the duration of recovery. • The power held during recovery interval is averaged and assumed to be constant. Fig. 4 : The power interval test protocol. After a warm-up period, the subject pedals at CP 4 for 2 minutes. Then the cyclist pedals at three different recovery power levels for different time intervals to recover energy. Following that, the subject performs a 3-min-all-out to burn all the remaining energy from AW C.
We have conducted experiments on a total of 17 subjects. All of the subjects are cyclists who cycle at least 3 times a week, and are used to high intensity workout sessions. Each subject was scheduled for 14 hour-long visits to our laboratory. Because of the complexities of such scheduling, only 6 of them were able to finish all of the tests. In order to be compliant with IRB policies we label our subjects by a number. The subjects who successfully finished all the experiments were Subs 6, 9, 11, 12, 14, 16. All tests were conducted in a laboratory setting on both Clemson and Furman University campuses. The tests were programmed on a RacerMate CompuTrainer [65] using Perfpro studio software application [66] . Figure 3 shows the experimental setup in our laboratory. On the first visit, the ramp test was conducted followed by a 3MT familiarization test. On four subsequent visits, a fresh 3MT was conducted. On the fifth visit, the interval test familiarization was conducted. On the subsequent visits, the interval tests at three different recovery powers (min power on Computrainer 80 Watts, 90% of power at GET (P GET ), and half way between P GET and CP ) and three durations (2 min, 6 min, 15 min) were conducted. There was at least 24 hours between two consecutive tests to ensure complete recovery. In the interval test protocol as shown in Figure 4 , after a warm up of 10 minutes, there is a 2 minute interval at CP 4 , the power at which all of the subject's AW C will be consumed in 4 minutes (CP 4 = AW C 240 ). The subject expends 50% of their AW C in the 2-min CP 4 interval and then recovers AW C at the mentioned recovery powers and durations. Following the recovery interval, the subject then performs a 3-min-all-out test to expend all of their remaining AW C. The amount of energy recovered in the recovery interval is then determined by subtracting AW C from the summation of areas above CP through the entire test.
The 3MT test results were used to plot remaining energy of each subject vs. their power during the test. Since the subjects are required and are constantly encouraged during the 3MT test to apply their maximum power, the mentioned plot can be used to find the parameter α of the model in Equation (5). As a sample, the experimental data from 3MT test for Sub 14 is presented in Figure 5 .
In order to find the recovery model parameters, the results of the interval tests are analyzed. As mentioned in Section III, it is assumed that the recovery model doesn't depend on time interval. Therefore, at each power level below CP , the adjusted power is calculated for the 3 time intervals and then averaged. By doing this we will have a plot of adjusted power vs actual applied power, consisting of 3 data points. Figure 6 represents the recovery model derived for Sub 14. The summary of modeling results for all of the subjects is presented in Table II .
In this section, the control problem formulation is discussed. First, the state space model of the system should be determined. The system has three states: i) traveled distance s ii) velocity of the bicycle v and iii) remaining energy of the cyclist w, and a single input which is rider's power u,
in which
where x(t) is the states vector, u(t) is the rider's power and the control input to the system. The function f is the nonlinear state-space model function, which relates the rate of change of the 3 states to the states and input. The first state function f 1 in Equation (6) is simply,
Newton's second law can be used to write the second state function f 2 which represents bicycle's acceleration:,
where u(t) is the applied power at the wheels, m t is the total mass of the bicycle and rider, C d is the drag coefficient, A is the frontal area, ρ is the density of air which is assumed to be constant and independent of the elevation, θ is the road slope which is positive for uphill and negative for downhill, and C R is the coefficient of rolling resistance of the road.
Assuming 100% efficiency for the bicycle powertrain, we can equate u(t) to the cyclist's power on the pedals. The advantage of using Equation (9) is that gear selection is not a factor in our formulation, which otherwise makes the optimization more complex. Note that the third state equation for time derivative of cyclist's remaining anaerobic energy is the previously represented Equation (3),
Now we are able to formulate a minimum-time optimal control problem to formalize the pacing strategy in a time-trial. Therefore, the objective function to be minimized is time,
The optimization constraints include the state-space model in Equation (6), limits on cyclist's power generation ability, the bicycle velocity, and energy shown below, velocity limits:
0 v(t) v max remaining energy limits: 0 w(t) AW C rider's power limit: 0 u(t) u max (t)
VI. ANALYTICAL APPROACH TO THE OPTIMAL CONTROL PROBLEM In this section we study the defined optimal control problem using the variational approach. According to Pontryagin's Minimum Principle (PMP), the necessary condition for the optimality of input u is that it minimizes the following Hamiltonian function,
where λ is the vector of co-state variables, L is the integrand in the cost function J in (11) , and f is the vector on the right hand side of the state equations with components represented in Equations (8), (9) , and (10). We also need to address the constraint of maximum power generation of the cyclist in Equation (15) . The maximum power u max is a function of state variable w as represented in (5) . Therefore, we observe an inequality constraint on a function of the control u and the third state w variables. We can rewrite the constraint in a standard form as,
In the presence of such inequality constraints on control and state variables, the Hamiltonian equation needs to be augmented as follows [67] ,
This will results in a switching Hamiltonian. However, the term µC will always be zero in Hamiltonian. The dynamics of the co-states should follow,
where the subscript x denotes the partial derivative with respect to the state variables. This is not the only switch in the Hamiltonian. Rider's remaining energy model in Equation (10) switches between fatigue and recovery conditions. It should be noted that in the recovery mode the constraint C is inactive, thus µ = 0. We only need to consider µ in the fatigue mode where the power level can reach its maximum and activate the constraint. A necessary condition for optimality is that the control input u, minimizes the Hamiltonian. One can set the partial derivative of H with respect to u equal to zero. However, in this case because H is affine in u, the derivative will be,
where subscript u denotes the derivative of H with respect to the input variable u. Also, the µ term only shows up when u CP because when u < CP the constraint C is not active. Since Equation (19) does not depend on u, the optimal solution will be of the bang-singular-bang form; that is the Hamiltonian is minimized at extreme values of u with the exception of potential singular arcs in between. The potential singular arcs can happen when the slope of the Hamiltonian (with respect to u) is zero.
When the Hamiltonian is affine in u, the sign of the line's slope indicates the optimal input value. As is shown in Equation (19), the slope depends on µ. When µ has a non-zero (positive) value, the constraint C in Equation (15) is active which means the optimal value for u is its maximum (u max ) regardless of the sign of the slope. Therefore, we only need to consider other cases where the constraint is inactive and µ = 0.
The system of equations that should be solved are the state-space Equations (6) and (18), which provide 6 equations combined. Let's open Equation (18),
where h(s) g (sin(θ(s)) + C R cos(θ(s))) (21)
We will assume the state w is free at final time, i.e. w(t f ) is free which requires that the corresponding co-state λ 3 be zero at the final time λ 3 (t f ) [67] . We do not force w to be zero and let the controller to choose the optimal trajectory at the end of the course. This means the boundary condition at the end of the course should be λ 3 = 0. On the other hand, according to Equation (20) , the rate of change of λ 3 depends on α and µ, both of which are non-negative values. Therefore, λ 3 is a negative value during the course which increases eventually to zero at the finish line. Moreover, now that we are considering cases where µ = 0 during the trial, λ 3 is a constant negative parameter.
Let's compare the slope of the Hamiltonian function in both fatigue and recovery modes. If the slope of the Hamiltonian line is positive, the minimum value of u minimizes the function. On the other hand, if the slope of the Hamiltonian line is negative, the maximum value of u minimizes the function. We can consider three cases,
• I: λ2 mv − λ 3 < 0 AND λ2 mv − aλ 3 < 0 The slope in both fatigue and recovery modes are negative which means maximum value of u in each case minimizes the Hamiltonian. Now, we should see which one results in a lower Hamiltonian. In the fatigue mode the maximum input value is u max , and in the recovery mode the maximum input value is CP ,
where the subscripts of H differentiates between fatigue and recovery modes because the equations of the two modes are slightly different. If we substitute the aforementioned values of u we cannot decisively say which Hamiltonian is smaller. We can write the optimal input as,
It should be noted that if the Hamiltonian is equal in both fatigue and recovery modes, we will have two optimal solutions. • II: λ2 mv − λ 3 > 0 AND λ2 mv − aλ 3 > 0 In this case, the input should take its minimum value in both cases, which will be CP and 0 for fatigue and recovery modes, respectively. The minimum Hamiltonian can be found from,
The optimal input value in this case will be:
Similar to the previous case, if the Hamiltonian is equal in both modes, we will have two optimal solutions. • III: λ2 mv − λ 3 > 0 AND λ2 mv − aλ 3 < 0 Since a has a constant value between 0 and 1, the slope of the Hamiltonian in recovery mode can be negative, while the slope in fatigue mode can be positive. The vice versa cannot be true. In this case, the input takes its minimum value in fatigue mode, and its maximum value in recovery mode. In both of these scenarios the input is CP , so we can decide,
So far the optimal control input can take values of 0, CP , and u max . This means, when needed, the cyclist should recover at zero power. A hypothetical recovery scenario can be when the cyclist is on a steep downhill. However, there is still another case that is not investigated, which is a possible singularity condition.
• Singular Arc: λ2 mv − λ 3 = 0 OR λ2 mv − aλ 3 = 0 Here we present the calculation for the case where the Hamiltonian's slope in the recovery or fatigue mode is zero. It can be shown that we get the exact same final result for both modes. Therefore, we proceed with the recovery mode. During the possible singular condition in recovery mode we have,
We can take the time derivative from both sides of the Equation (27) . From Equation (20) we know that the time derivative of λ 3 is zero,
We can substituteλ 2 andv with Equations (20) and (9), respectively:
Simplifying the equation above by using (27) yields,
Since input u does not appear in Equation (30), we need to take the time derivative one more time,
We can simplify the above equation using Equation (20) . Then,
In Equation (32) the only parameter that can be zero isv. Therefore, during the singular interval the velocity must be a constant. The corresponding input power is obtained using Equation (9),
which varies with the road grade.
Considering all of the cases discussed above, the optimal power trajectory can only take values from the vector below,
Thus we have succeeded in identifying the optimal modes in a time-trial strategy. That is, the cyclist either outputs maximal power, comes to standstill to rest, rides at CP, or else must pedal at a constant speed which means variable power if the road grade is varying. While this insight is useful, it is difficult to analytically determine when the switching between these modes happens. Note that the three state equations along with the three co-state equations form a two point boundary value problem that is generally difficult to solve analytically. The switching dynamics of the problem at hand creates additional challenges. Therefore, we resort to numerical solution of the optimal control problem via dynamic programming. The insights obtained from the analytical treatment could help by limiting the input space when iterating on DP solutions.
In our DP implementation, at every step in distance, the state and input variables are discretized between their minimum and maximum. First, the state space model in Equation (6) is discretized. Since the final destination is known, we pick distance as the independent variable and discretize it at regular intervals ∆s = 100 m. With a zero-order hold on inputs in between sampling intervals we obtain the following discretized state-space equations:
The cost function in Equation (11) is rewritten with position as the independent variable and discretized as follows,
According to the Bellman's principle of optimality [68] , when a system is on its optimal path from an initial state to a final state, regardless of any past decision or state, it should follow an optimal policy for the remainder of the route. Therefore, in dynamic programming, to find the optimal state trajectory, one can begin from the final state and move backward and calculate the optimal cost-to-go from any state to the final step. Figure 7 demonstrates the backward dynamic programming method using the Bellman's principle of optimality. The optimal costs from all of the nodes at s i+1 to the final state at s N are stored as J * i+1,N . Then, cost to go from every node at s i to all of the nodes at s i+1 is calculated and the optimal value among them is stored at the specified node at s i . This process should be continued backward to the beginning of the route. Subsequently, in a forward DP sweep, the optimal action at each discretized state node will be known. Let's denote the optimal cost-to-go from specific velocity and energy states at step s i+1 to s N by J * i+1,N . Then, the optimal trajectory from step s i to s N will be,
where
Initially, the cyclist starts from zero velocity and his/her remaining energy is initialized at AW C.
The DP formulation was then coded in MATLAB. At every step in distance, velocity v is quantized into 32 nodes from 0.1 to 16 m s . Anaerobic energy w is quantized into 100 nodes from 0 to the Sub 14th AW C (7841J). Input power u is also quantized by 298 nodes from 0 to 594 Watts. When at the distance s i from the initial position, v i and w i states move to v i+1 and w i+1 by applying an input u i . The resultant states at step s i+1 will not necessarily attain the quantized values of states v and w. Therefore, we need to interpolate between every node in the velocity-energy plane. Under this condition we select the closest available node to the calculated states at step s i+1 . 1 This was to avoid a time consuming process of re-quantizing the states at every step. The DP simulation was ran on a laptop with a 1.3GHz Intel core i7 CPU, and 16GB of RAM. The run time for our longest cycling route (22km) was 35.8 minutes.
In this set of simulations, we have chosen three courses. First, a flat road is considered to observe the optimal strategy when road grade is zero. Second, the cycling course of the 2019 Duathlon National Championship in Greenville, South Carolina, is simulated which has mild variable road grade. Lastly, a difficult cycling route in Leysin, Switzerland, is picked. The elevation map for both real courses are downloaded from Strava [69] and shown in Figure 8 . Figure 9 shows the result of the optimal pacing simulations with Sub 14th data. On the flat road, the optimal pace is to go all out and then coast at CP . On the Greenville course, the rider should go all-out in the beginnig but reserve a little energy for the upcoming hills one the course. During the rest of the course, either the maximum power constraint is activated, the subject is pedaling at constant speed during singular intervals, or power is at CP . These values are in accordance with the optimal control actions in Equation (34) . The most interesting outcome was from simulating the Switzerland course. We can observe bang-bang control between maximum power and zero. The optimal strategy is to perform several all-out intervals with enough rest in between to recover the entire AW C. These all-out intervals are needed because in several sections of the course, the road grade is too high for Sub 14 to climb with any power lower than her absolute maximum. Therefore, there are multiple situations where the cyclist has to stop pedaling to recover her anerobic work capacity. The cyclist stops pedaling for 400 sec during which the cyclist can recover all of her anaerobic energy. This means riding at CP does not provide enough power to overcome the steep hills in this route. In fact, this course is not a standard course for a time trial since usually the average road grade for a time trial race should be between 1% to 6% [70] . The road grade in this course averages at 11% and goes up to 20% at several intervals. Hence, this simulation just has an educational benefit of validating our analytical approach in Section VI.
The above simulations indicate that, the general strategy for an individual time-trial effort is to start with maximum power and then try to pedal in a close proximity of CP and at a constant velocity, unless there are hills with high grade. In the cycling community, it is commonly discussed that an effective time trial strategy includes riding at a relatively constant power [70] . Our simulations with zero and moderate road grades prove that this is actually the mathematically optimal strategy. That said, the key factor is selecting the right power level for the course, which our optimal control solution provides.
We had the opportunity to simulate the Greenville Duathlon course on our Computrainer. 2 We ran this experiment with Sub 14 by requesting the subject to ride with her own strategy as the baseline test. During the test, she was receiving velocity, 2 One of the shortcomings of utilizing Computrainer is that it does not model drag force during a ride. Therefore, the drag term in Equation (36) should be removed. Also, when the subject is riding on a downhill, Computrainer cannot accelerate the bicycle as would happen on a real road. Therefore, besides removing the negative resistance force due to downhill section of the course, the calibrated resistance force of 15.5 N is added to Equation (36) instead of the drag force. This force is imposed constantly by the computrainer and is set during the calibration process. distance, power, elevation, and heart rate data in real-time via the interface depicted in Figure 3 . The subject managed to finish the course in 34 min 8 sec. However, our DP simulation suggests that if she followed the optimal power provided by the DP, she would have finished the course in 27 min 41 sec. An important observation from this test is that the subject's self strategy lacks consistency of applied power. The standard deviation of power data, shown in Figure 10 , for the simulation is 7.2 Watts lower than the baseline test. Also, the energy curve in Figure 10 suggests that the cyclist was pedaling below CP during the second half of the course. Therefore, she recovered all of her anaerobic energy at the end which is not optimal; one should finish a time-trial completely depleted as the optimal w trajectory in Figure 10 shows.
After completing the simulations in Section VIII-A and to test versus a commonly recognized baseline, we simulated our Sub 14 in a virtual stage of Tour de France competition. Tour de France is commonly described as the most prestigious and difficult cycling race in the world [71] . The best cyclists around the world attend this 23-day event. The competition consists of 21 stages. In 2019, the 13th stage of Tour de France was an individual time trial. In the section, we simulate the stage 13 of the competition to find out what the optimal strategy would be for our Sub 14 if she participated in Tour de France. Unfortunately, Tour de France is a male-only competition, and Sub 14 is a semi-professional female cyclist. Therefore, we are not comparing Sub 14 to the Tour de France athletes. We want to gain insight into the general strategy over the course.
As discussed earlier in Section VIII-A, we showed that the common strategy of maintaining a relatively constant power during a time trial effort is inline with our optimal control results. However, the optimal strategy can slightly vary depending on the elevation profile of a course. On the elevation plot of the stage 13 time trial competition presented in Figure  11 , three downhill sections are highlighted in red boxes. The optimal power trajectory in Figure 11 shows that the cyclist should stop pedaling in the downhill sections to recover some energy. Then, at the end of the downhill, the maximum power constraint is activated and the cyclist is supposed to perform a short all-out interval. This is in contrast with the constant power strategy commonly used by cyclists.
According to the 2019 Tour de France ranking at stage 13, our Sub 14 would finish last among the competitors. The stage winner finished the course in 35 minutes, while our simulation suggests Sub 14 could finish in 44 minutes. She is 2 min and 5 seconds behind the last person who finished the race in 2019. It is important to note that this time can be achieved in ideal condition where it is assumed that the cyclist could exactly follow the optimal power trajectory.
In this study, we proposed a set of mathematical models for cyclists' fatigue and recovery during anaerobic exercise. Three experimental protocols were designed to calculate the model parameters which were done using a Computrainer.The modeling results for 6 subjects are presented. Then, these models are used to formulate the pacing strategy in a time-trial case as an optimal control problem. The analytical investigation of the problem suggests a bang-singular-bang solution. Since the analytical solution is too complex to be explicitly calculated, the dynamic programming method is employed for numerical solution of optimal pacing. We simulated three courses to make sure that DP results are in accordance with our analytical solution. Another experimental protocol is designed to simulate the Greenville Duathlon course on the Computrainer. Sub 14 was instructed to ride on the Computrainer with their self strategy. Then we used the test results as a baseline to compare with our optimal simulation. The simulation suggested that Sub 14 could finish the course 6 min and 27 sec sooner if the optimal power was followed. The comparison shows that consistency in power generation is the key difference between the subject's self and the optimal strategies. As a motivation for future work, we simulated the 2019 Tour de France stage 13 which was an individual time trial competition. Although our subject finished last in the race, the results provide valuable insights about the strategy for the time-trial stage. For example, sometimes it is optimal to stop pedaling in downhills to recover some energy to overcome an upcoming hill. This strategy can be counterintuitive.
There is still a lot of room for improvement. The most important of all, we have set up a test protocol for providing optimal power to our subjects in real time, and observe their performance improvement. At the time of writing this paper, unfortunately we are forced to pause our testing because of the spread of COVID-19, and may have to repeat the calibration tests or with new subjects in a few months time. Eventually, we would like to perform a outdoor tests on a real course. For achieving this goal, we need to instrument a bicycle with necessary sensors and equipment which is the future plan of our group. 