电动无人方程式赛车时间最优控制方法研究

Research on time-optimal control method for an autonomous electric fomula race car

摘要

针对电动无人方程式赛车高速循迹场景，提出一种最小化单圈完赛时间的最优控制方法.基于赛车的纵向动力学和动力系统模型，构建时间域内的最短时间最优控制问题，结合非光滑分析与庞特里亚金极小值原理推导得到最优控制问题的解析控制策略.为获取解析控制策略中的关键协态变量，将时间域最优控制问题转化为位置域最优控制问题，并进一步将其转化为二阶锥形式的凸优化问题.采用基于MATLAB的YALMIP和MOSEK求解工具对该凸优化问题进行建模和数值求解，获取关键协态变量的最优轨迹.模拟计算结果表明，结合凸优化给出的关键协态变量，解析控制策略可给出准确的最优驱动和制动操作.

Abstract

Focusing on high-speed tracking scenarios of autonomous electric formula race cars， this paper proposes an optimal control method for minimizing the completion time of one lap. Firstly， a minimum-time optimal control problem is constructed in the time domain based on the longitudinal dynamics and driving system model of a race car. By combining nonsmooth analysis with Pontryagin’s minimum principle， an analytical optimal control strategy is derived. Secondly， to obtain the key costates of the analytical optimal control strategy， the optimal control problem is transformed from the time domain to the space domain， and further converted into a convex optimization problem in the second-order cone form. Finally， the convex optimization problem is modeled and numerically solved using MATLAB’s YALMIP and MOSEK solver tools， and therefore， the optimal trajectories of key costates are obtained. Simulation results show that， with the key costates obtained from convex optimization， the analytical optimal control strategy can provide accurate optimal driving and braking operations.

Graphical abstract

关键词

最优控制 / 非光滑分析 / 庞特里亚金极小值原理 / 凸优化 / 无人方程式赛车

Key words

optimal control / nonsmooth analysis / Pontryagin’s minimum principle / convex optimization / autonomous formula race car

无人驾驶是汽车未来发展的关键方向，汽车全面智能化的发展浪潮推动了众多企业进军这些领域^［1］.为了培养无人驾驶与智能网联汽车行业的综合性人才，中国汽车工程学会在大学生电动方程式大赛的基础上，进一步推出了中国大学生无人驾驶方程式大赛（formula student autonomous China， FSAC）.高速循迹赛是大赛项目的重要组成部分，无人方程式赛车需要在满足各项要求的前提下以最短时间完成预定赛道.对于纯电动赛车而言，电耗是限制车辆圈速表现的关键因素，因此如何在车辆动力学以及电功率的限制下进行车辆圈速的最优化成为各电动车队关注的重点问题^［2］.

国内外研究人员对于纯电动赛车的最快圈速优化问题进行了诸多研究.Sciarretta等^［3］全面阐述了车速优化的典型最优控制问题形式，描述了这类问题的通用数学模型，以及常用动态规划法、直接法与间接法.Dong等^［4-5］考虑道路形状特征和限速等信息，使用动态规划法在一条长为8 km的仿真道路上求得了纯电动车辆的最优车速轨迹，但计算时间成本较高.为降低计算时间成本，其通过迭代的形式优化动态规划算法的流程，每次迭代时可自适应地缩小变量搜索范围.Kim等^［6］研究了基于动态规划的燃料电池动力系统赛车速度曲线和功率分配的协同优化方法，并应用庞特里亚金极小值原理减少控制量，降低了动态规划的计算成本.Herrmann等^［7］基于无人驾驶电动赛车的传动系功率损失和电机热模型建立最快圈速优化的非线性规划问题，并通过数值求解得到赛车在热力学限制下的最优比赛策略.Padilla等^［8］对全局车速优化问题进行了重新建模，引入速度相关的附加项并进行合适的离散化，给出了一种求得全局最优解的序列二次规划算法.Hamednia等^［9］提出了一种高效的分层车速优化算法，上层在线求解以行驶时间和车速为状态变量的非线性规划问题，下层离线解决挡位优化的问题，在模型预测控制框架下应用序列二次规划算法求解了最优车速.Li等^［10］总结了凸优化方法在各类电动车零部件尺寸优化和能量管理问题中的应用，讨论了凸优化方法在电动车设计与控制中的应用前景和发展趋势.Broere等^［11］考虑赛车单个轮胎的抓地力限制，结合车辆动力学与电机、电池的凸模型，提出了一种计算全轮驱动纯电动赛车最快圈速控制策略的凸优化数值求解方法，可高效求解全局最优的驱动电机扭矩.Borsboom等^［12］基于不同类型的赛车变速箱与电机模型建立了电动赛车最快圈速的凸优化问题，并通过二次规划算法数值求解了全局最优的电机扭矩，证明了在某一预定赛道上CVT变速箱的优势.Van Kampen等^［13］考虑车辆动力总成部件的热力学限制，提出了一种针对耐力赛车的双层混合整数凸优化框架，上层为单圈行驶时间优化，下层为行驶距离优化.De Castro等^［14］基于非线性四轮驱动电动车模型，通过凸优化的方法处理轮胎力与功率约束，设计了基于最短时间问题的路径跟踪控制器.Salazar等^［15-17］考虑赛规要求与抓地力限制，结合车辆动力学与功率流模型离线推导了单圈最短时间问题最优控制策略的解析解，在此基础上设计了基于三个PID控制器的反馈控制策略，通过控制关键变量的值实时跟踪最优控制策略，并最终提出了一种基于两级模型预测控制的实时控制算法用于定期更新最优控制策略.Cartignij等^［18］在Salazar的基础上考虑多挡变速器与换挡策略对最短时间问题的影响.Duhr等^［19］进一步研究了车辆电池能耗约束对于解析形式最优控制策略的影响.

综上可以看出，赛车最快圈速优化问题的研究总体上可分为直接的数值求解与间接的解析求解两种思路，数值求解可基于合理的问题建模使用高效的求解器在短时间内获得关键变量的解，解析求解可利用最优化理论推导出在限定条件下精确的解析形式最优控制策略.本文结合解析解的理论精确性与数值解的求解高效性，提出了最短时间问题最优控制策略的推导与求解方法.本文主要包括三方面内容.第一部分内容为在时间域内构建时间最短的最优控制问题，结合非光滑分析^［20］与庞特里亚金极小值原理^［21］解析推导出与关键协态变量相关的最优控制策略函数；第二部分内容为在位置域内将时间最短的最优控制问题重新建模为二阶锥形式^［22］的凸优化问题^［23］；最后一部分内容为通过MATLAB的高效求解工具进行数值求解得到关键协态变量的变化轨迹，进一步求得赛车在赛道每个位置的最优控制量解析解结果，并通过解析解的模拟实施结果与高效求解器数值结果对比验证了该最优控制策略的合理性.

1 单圈时间最优控制策略函数

1.1 车辆纵向动力学

本文以某无人方程式赛车作为研究对象，车辆的纵向动力学主要考虑其动能与位置变化，动能E_{\mathrm{k}\mathrm{i}\mathrm{n}}(t)变化在时间域内表示为：

\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}(t)=P_{\mathrm{p}}(t)-P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(t)-P_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}(t))（1）

式中：P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(t)为制动功率；P_{\mathrm{p}}(t)为车辆的轮上驱动功率，其将在1.2节中被进一步定义；P_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}(t))为包括空气阻力与滚动摩擦阻力在内的总阻力功率，其定义见式（2）：

P_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}(t))=(\mathrm{0.5}A\rho C_{\mathrm{d}}\frac{\mathrm{2}}{m}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}+mg{\mu }_{\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{l}})\sqrt[]{\frac{\mathrm{2}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{m}}（2）

式中：A为车辆迎风面积；\rho为空气密度；C_{\mathrm{d}}为风阻系数；m为车重；g为重力加速度；{\mu }_{\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{l}}为滚动摩擦系数.本文参考赛道的道路坡度角恒为零.

基于上述参数可得车辆纵向位置s(t)的变化为：

\frac{\mathrm{d}}{\mathrm{d}t}s(t)=\sqrt[]{\frac{\mathrm{2}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{m}}（3）

1.2 车辆动力系统

一辆无人驾驶纯电动赛车的基础动力系统功率流如图1所示.

图1中：P_{\mathrm{b}\mathrm{a}\mathrm{t}}为电池输出功率；P_{\mathrm{b}\_\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}为电池内部损失的功率；P_{\mathrm{p}}、P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}、P_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}的定义与前文相同；P_{\mathrm{i}}为电池总功率，其与电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}的关系见式（4）：

P_{\mathrm{i}}(t)={\alpha }_{\mathrm{b}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}{(t)}^{\mathrm{2}}+P_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)（4）

式中：{\alpha }_{\mathrm{b}}\ge \mathrm{0}为电池效率参数.此关系式中考虑了P_{\mathrm{b}\_\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}的影响.通过一阶RC电路模型可推导得：

P_{\mathrm{i}}=P_{\mathrm{b}\mathrm{a}\mathrm{t}}+{\left(\frac{U_{\mathrm{o}\mathrm{c}}-\sqrt[]{{U_{\mathrm{o}\mathrm{c}}}^{\mathrm{2}}-\mathrm{4}R{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}}}{\mathrm{2}R}\right)}^{\mathrm{2}}R（5）

式中：U_{\mathrm{o}\mathrm{c}}为电池开路电压；R为电池内阻.进一步结合\sqrt[]{{U_{\mathrm{o}\mathrm{c}}}^{\mathrm{2}}-\mathrm{4}R{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}}的二阶泰勒展开式^［24］可得到：

P_{\mathrm{i}}=\frac{R}{{U_{\mathrm{o}\mathrm{c}}}^{\mathrm{2}}}{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{2}}+P_{\mathrm{b}\mathrm{a}\mathrm{t}}（6）

即电池效率参数{\alpha }_{\mathrm{b}}等于R/{U_{\mathrm{o}\mathrm{c}}}^{\mathrm{2}}.

电池总能耗在时间域的变化率为电池总功率：

\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{b}}(t)=P_{\mathrm{i}}(t)（7）

考虑电机及传动系的效率，定义车辆轮上驱动功率P_{\mathrm{p}}(t)与P_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)的关系为：

P_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)={\alpha }_{\mathrm{p}}{P}_{\mathrm{p}}{(t)}^{\mathrm{2}}+P_{\mathrm{p}}(t)（8）

式中：{\alpha }_{\mathrm{p}}\ge \mathrm{0}为由数据拟合得到的电机及传动系效率参数^［25］.拟合效果如图2所示.

1.3 变量约束

动能约束来自赛道不同路段的曲率，表示为：

\mathrm{0}\le E_{\mathrm{k}\mathrm{i}\mathrm{n}}(t)\le E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s(t))（9）

式中：E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s(t))为车辆在位置s处的最大动能.

对于单圈优化问题，期望车辆在单圈比赛末时刻T的速度与起步时相同，用动能表示为：

E_{\mathrm{k}\mathrm{i}\mathrm{n}}(T)=E_{\mathrm{k}\mathrm{i}\mathrm{n}}(\mathrm{0})（10）

赛车完成单圈的电池总能耗约束取决于不同的比赛策略，表示为：

E_{\mathrm{b}}(T)\le E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}（11）

式中：E_{\mathrm{b}}(T)是赛车完成单圈赛道的电池总能耗；E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}是允许的单圈最大电池能耗，可根据不同的比赛策略选定不同的值.

最后，控制量电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)与制动功率P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(t)的约束分别为：

P_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)\in [\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}]（12）

P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(t)\in [\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}]（13）

式中：P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}和P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}分别是最大电池输出功率与最大制动功率.

1.4 最短单圈时间最优控制问题

1.4.1 目标函数与状态方程

时间最优控制问题的目标是使单圈时间T最短.定义状态变量为车辆动能、电池总能耗、车辆已行驶路程，\mathit{x}=(E_{\mathrm{k}\mathrm{i}\mathrm{n}},E_{\mathrm{b}}{,s)}^{\mathrm{T}}.控制量为电池输出功率、制动功率，\mathit{u}=(P_{\mathrm{b}\mathrm{a}\mathrm{t}},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}{)}^{\mathrm{T}}.完整的最优控制问题表述如下：

\mathrm{m}\mathrm{i}\mathrm{n}{\displaystyle \underset{\mathrm{0}}{\overset{t}{\int }}}g(\mathit{x},\mathit{u})dt+h(x(T))

\left\{\begin{array}{l}{\dot{E}}_{\mathrm{k}\mathrm{i}\mathrm{n}}=P_{\mathrm{p}}-P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}-P_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}})\\ {\dot{E}}_{\mathrm{b}}=P_{\mathrm{i}}\\ \dot{s}=\sqrt[]{\frac{\mathrm{2}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{m}}\end{array}\right.

\left\{\begin{array}{l}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)\in [\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}]\\ P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(t)\in [\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}]\end{array}\right.

\left\{\begin{array}{l}\mathrm{0}\le E_{\mathrm{k}\mathrm{i}\mathrm{n}}(t)\le E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s(t))\\ E_{\mathrm{k}\mathrm{i}\mathrm{n}}(T)=E_{\mathrm{k}\mathrm{i}\mathrm{n}}(\mathrm{0})\\ E_{\mathrm{b}}(T)\le E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}\\ s(T)=S\end{array}\right.（14）

式中：S为单圈赛道的总长度；T为赛车完成单圈赛道所需的总时间；其中P_{\mathrm{p}}与P_{\mathrm{i}}的定义与前文均相同，如式（15）所示：

\left\{\begin{array}{l}{P}_{\mathrm{p}}=\frac{\sqrt[]{\mathrm{4}{\alpha }_{\mathrm{p}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}+\mathrm{1}}}{\mathrm{2}{\alpha }_{\mathrm{p}}}-\frac{\mathrm{1}}{\mathrm{2}{\alpha }_{\mathrm{p}}}\\ P_{\mathrm{i}}={\alpha }_{\mathrm{b}}{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{2}}+P_{\mathrm{b}\mathrm{a}\mathrm{t}}\end{array}\right.（15）

目标函数中的g(x,u)为过程代价函数，其定义为：

g(x,u)=\mathrm{1}+{\psi }_{\left[-\mathrm{\infty },\mathrm{0}\right]}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}-E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s))+

{\psi }_{\left[\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\right]}(P_{\mathrm{b}\mathrm{a}\mathrm{t}})+{\psi }_{\left[\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}\right]}(P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}})（16）

目标函数中的h(x(T))为终端代价函数，其定义为：

h(x(T))={\psi }_{\left\{E_{\mathrm{k}\mathrm{i}\mathrm{n}}(\mathrm{0})\right\}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}(T))+{\psi }_{[-\mathrm{\infty },E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}]}(E_{\mathrm{b}}(T))+

{\psi }_{\left\{S\right\}}(s(T))（17）

式（16）和式（17）中的{\psi }_{\mathrm{\xi }}(x)函数为非光滑的指示函数，其定义为：

{\psi }_{\xi }(x)=\left\{\begin{array}{cc}\mathrm{0},& x\in \xi \\ \mathrm{\infty },& x\notin \xi \end{array}\right.（18）

式中：\xi为指示函数的定义域.引入指示函数可在最优控制问题中简化不等式约束的处理.指示函数的定义域\xi为凸集，其在x\in \xi处的偏微分是集合\xi的法锥，定义为：

\frac{\partial }{\partial x}{\psi }_{\xi }(x)=N_{\xi }(x)=

\{y\in R^n|{\mathit{y}}^{\mathrm{T}}(\overline{x}-x)\le \mathrm{0},x\in \xi ,\overline{x}\in \xi \}（19）

式中：N_{\xi }(x)即为集合\xi的法锥；{\mathit{y}}^{\mathrm{T}}表示向量 y 的转置.已知一个凸集[a,b]\in R，则以该集合为定义域的指示函数及法锥如图3所示.

1.4.2 最优控制策略函数推导

基于前文所述的目标函数与状态方程，可得哈密顿函数为：

H(x,u,\lambda )=\mathrm{1}+{\lambda }_{\mathrm{1}}\left(\frac{\sqrt[]{\mathrm{4}{\alpha }_{\mathrm{p}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}+\mathrm{1}}}{\mathrm{2}{\alpha }_{\mathrm{p}}}-\frac{\mathrm{1}}{\mathrm{2}{\alpha }_{\mathrm{p}}}-\right.

\left.P_{brake}-P_{drag}(E_{kin}{)}^{\begin{array}{l}\\ \\ \end{array}}\right)+{\lambda }_{\mathrm{2}}({\alpha }_{\mathrm{b}}{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{2}}+P_{\mathrm{b}\mathrm{a}\mathrm{t}})+

{\lambda }_{\mathrm{3}}\sqrt[]{\frac{\mathrm{2}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{m}}+{\psi }_{\left[-\mathrm{\infty },\mathrm{0}\right]}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}-E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s))+

{\psi }_{\left[\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\right]}(P_{\mathrm{b}\mathrm{a}\mathrm{t}})+{\psi }_{\left[\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}\right]}(P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}})（20）

式中：{\lambda }_{\mathrm{1}}~{\lambda }_{\mathrm{3}}为分别对应于状态量E_{\mathrm{k}\mathrm{i}\mathrm{n}}、E_{\mathrm{b}}、s的协态变量，其代表了改变各状态量所需的最优代价.其余变量的定义与前文相同.

各协态变量在时间域的微分方程可根据下式求得：

\stackrel{•}{{\lambda }_i}=-\frac{\partial }{\partial x_i}H(x^{\mathrm{*}},u^{\mathrm{*}},\lambda )\begin{array}{cc}& i=\mathrm{1,2},\mathrm{3}\end{array}（21）

可得：

\stackrel{•}{{\lambda }_{\mathrm{1}}}\in \frac{\partial P_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}}{\partial E_{\mathrm{k}\mathrm{i}\mathrm{n}}}({E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{}}^{\mathrm{*}}){\lambda }_{\mathrm{1}}-\frac{{\lambda }_{\mathrm{3}}}{\sqrt[]{\mathrm{2}m{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}}}}-

N_{\left[\mathrm{0},E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s)\right]}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}})\stackrel{•}{{\lambda }_{\mathrm{2}}}=\mathrm{0}

\stackrel{•}{{\lambda }_{\mathrm{3}}}\in -N_{\left[{E_{\mathrm{k}\mathrm{i}\mathrm{n}}}^{\mathrm{*}},+\mathrm{\infty }\right]}(E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s_{}^{\mathrm{*}}))\frac{\partial E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s_{}^{\mathrm{*}})}{\partial s}（22）

各协态变量在单圈末时刻T时的值可根据下式求得：

{\lambda }_i(T)=\frac{\partial }{\partial x_i}h(x^{\mathrm{*}}(T))\begin{array}{cc}& i=\mathrm{1,2},\mathrm{3}\end{array}（23）

可得：

\stackrel{}{{\lambda }_{\mathrm{1}}(T)}\in N_{\left\{{E_{\mathrm{k}\mathrm{i}\mathrm{n}}}^{\mathrm{*}}(\mathrm{0})\right\}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}}(T))

\stackrel{}{{\lambda }_{\mathrm{2}}}(T)\in N_{[-\mathrm{\infty },E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}]}(E_{\mathrm{b}}^{\mathrm{*}}(T))

\stackrel{}{{\lambda }_{\mathrm{3}}(T)}\in N_{\{S\}}(s_{}^{\mathrm{*}}(T))（24）

结合式（22）与式（24），可得{\lambda }_{\mathrm{2}}全程为常量且非负，当且仅当满足{E_{\mathrm{b}}^{}}^{\mathrm{*}}(T)=E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}时可取非零值.{\lambda }_{\mathrm{1}}与{\lambda }_{\mathrm{3}}均为时变量.

根据庞特里亚金极小值原理，可通过哈密顿函数推导最优控制策略：

\mathrm{0}\in \frac{\partial H(x^{\mathrm{*}},u^{\mathrm{*}},\lambda )}{\partial P_{\mathrm{b}\mathrm{a}\mathrm{t}}}=\frac{\mathrm{1}}{\sqrt[]{\mathrm{4}{\alpha }_{\mathrm{p}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}+\mathrm{1}}}{\lambda }_{\mathrm{1}}+

{\lambda }_{\mathrm{2}}(\mathrm{2}{\alpha }_{\mathrm{b}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}+\mathrm{1})+N_{[\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t}.\mathrm{m}\mathrm{a}\mathrm{x}}]}({P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}})（25）

\mathrm{0}\in \frac{\partial H(x^{\mathrm{*}},u^{\mathrm{*}},\lambda )}{\partial P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}}=-{\lambda }_{\mathrm{1}}+N_{[\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}.\mathrm{m}\mathrm{a}\mathrm{x}}]}(P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}})（26）

由式（26）可以进一步推导出：如果{\lambda }_{\mathrm{1}}为正数，则法锥N_{[\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}.\mathrm{m}\mathrm{a}\mathrm{x}}]}(P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}})一定为正数，因此P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}只能取得定义域区间的右端点，即P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}=P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}.同理，{\lambda }_{\mathrm{1}}为负数时，P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}只能取得定义域区间的左端点，即P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}=\mathrm{0}.{\lambda }_{\mathrm{1}}为0时，P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}\in [\mathrm{0},P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}].则最优制动功率P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}与{\lambda }_{\mathrm{1}}的完整关系可表示为：

（27）

由式（25）可通过类似的思路推导得到最优电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}与{\lambda }_{\mathrm{1}}的部分关系为：

P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})=\left\{\begin{array}{cc}\begin{array}{l}{P}_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\\ \mathrm{0}\end{array}& \begin{array}{l},{\lambda }_{\mathrm{1}}\le {\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\\ ,{\lambda }_{\mathrm{1}}\ge {\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}}\end{array}\end{array}\right.（28）

式中参数定义如下：

\left\{\begin{array}{l}{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}=-{\lambda }_{\mathrm{2}}(\mathrm{2}{\alpha }_{\mathrm{b}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}+\mathrm{1})\sqrt[]{\mathrm{4}{\alpha }_{\mathrm{p}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}+\mathrm{1}}\\ {\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}}=-{\lambda }_{\mathrm{2}}\end{array}\right.（29）

参数{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}与{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}}分别是最优电池输出功率{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}}取得最大值与最小值的临界{\lambda }_{\mathrm{1}}值，这两个值均与{\lambda }_{\mathrm{2}}直接相关.

在{\lambda }_{\mathrm{1}}\in ({\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}},{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}})的区间内，{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}}需通过另外的方法推导.{\lambda }_{\mathrm{1}}\in ({\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}},{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}})时，{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}}会取到上下限之间的值，即{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}}\in (\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}})，其对应的法锥N_{[\mathrm{0},P_{\mathrm{b}\mathrm{a}\mathrm{t}.\mathrm{m}\mathrm{a}\mathrm{x}}]}({P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}})=\mathrm{0}，则根据式（25）可推导出：

{\lambda }_{\mathrm{1}}=-{\lambda }_{\mathrm{2}}(\mathrm{2}{\alpha }_{\mathrm{b}}{P}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}+\mathrm{1})\sqrt[]{\mathrm{4}{\alpha }_{\mathrm{p}}{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}}+\mathrm{1}}（30）

进一步可得：

\mathrm{0}=\mathrm{16}{{\alpha }_{\mathrm{b}}}^{\mathrm{2}}{\alpha }_{\mathrm{p}}(P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}{)}^{\mathrm{3}}+(\mathrm{16}{\alpha }_{\mathrm{b}}{\alpha }_{\mathrm{p}}+\mathrm{4}{{\alpha }_{\mathrm{b}}}^{\mathrm{2}})(P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}{)}^{\mathrm{2}}+

(\mathrm{4}{\alpha }_{\mathrm{b}}+\mathrm{4}{\alpha }_{\mathrm{p}})P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}+\mathrm{1}-\frac{{{\lambda }_{\mathrm{1}}}^{\mathrm{2}}}{{{\lambda }_{\mathrm{2}}}^{\mathrm{2}}}（31）

式（31）解析地表达了{\lambda }_{\mathrm{1}}\in ({\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}},{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}})时对应的P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}，则最优电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}与{\lambda }_{\mathrm{1}}的完整关系为：

P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})=\left\{\begin{array}{cc}\begin{array}{l}{P}_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\\ \mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}(\mathrm{31})\\ \mathrm{0}\end{array}& \begin{array}{l},{\lambda }_{\mathrm{1}}\le {\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\\ ,{\lambda }_{\mathrm{1}}\in ({\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}},{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}})\\ ,{\lambda }_{\mathrm{1}}\ge {\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}}\end{array}\end{array}\right.（32）

式中：root（31）表示式（31）的根.最终所得的解析形式最优控制策略曲线{P_{\mathrm{b}\mathrm{a}\mathrm{t}}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})与{P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})与图4类似.

由图4可见最优控制策略为协态变量{\lambda }_{\mathrm{1}}的分段函数，且此函数的转折点与{\lambda }_{\mathrm{2}}直接相关.在式（20）中，{\lambda }_{\mathrm{1}}与控制量P_{\mathrm{b}\mathrm{a}\mathrm{t}}、P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}直接相乘，{\lambda }_{\mathrm{2}}与控制量P_{\mathrm{b}\mathrm{a}\mathrm{t}}直接相乘，因此{\lambda }_{\mathrm{1}}与{\lambda }_{\mathrm{2}}是影响最优控制策略的关键协态变量，图4也反映了这一点.{\lambda }_{\mathrm{3}}虽为时变量，但其在式（20）中并未与任何一个控制量直接相乘，因此最优控制策略并不与{\lambda }_{\mathrm{3}}直接相关，图4的函数曲线也并未体现{\lambda }_{\mathrm{3}}的作用.

2 位置域凸优化问题

前文1.4.2节已在时间域推导出与关键协态变量{\lambda }_{\mathrm{1}}、{\lambda }_{\mathrm{2}}直接相关的最优控制策略函数，为获得实际可参考的最优控制量，需求得{\lambda }_{\mathrm{1}}与{\lambda }_{\mathrm{2}}的最优值.凸优化问题是一种理论相当成熟的数学优化问题，通过高效的求解器，此类问题能快速地找到满足精度要求的全局最优解且具有很高的鲁棒性^［26］.因此，对于各类约束条件已定的特定场景单圈比赛，可将最优控制问题转化为特定形式的凸优化问题进行数值求解，从而获得{\lambda }_{\mathrm{1}}与{\lambda }_{\mathrm{2}}的单圈全程变化轨迹数值解.本文的优化目标是最短化单圈比赛的完成时间，总时间T为未知量，若在时间域内进行数值求解则会面临离散化的合理性问题，因此本文的凸优化数值求解在位置域内进行.

2.1 位置域优化问题构建

将前文1.4.1节的时间域优化问题转化到位置域，对应可得优化问题为：

\mathrm{m}\mathrm{i}\mathrm{n}{\displaystyle \underset{\mathrm{0}}{\overset{S}{\int }}}\frac{\mathrm{d}t}{\mathrm{d}s}(s)\mathrm{d}s

\left\{\begin{array}{l}\frac{\mathrm{d}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}\left(s\right)}{\mathrm{d}s}=F_{\mathrm{p}}(s)-F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(s)-F_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s))\\ \frac{\mathrm{d}{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}\left(s\right)}{\mathrm{d}s}=F_{\mathrm{i}}(s)\end{array}\right.

\left\{\begin{array}{l}\mathrm{0}\le F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)\le P_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}\frac{\mathrm{d}t}{\mathrm{d}s}(s)\\ \mathrm{0}\le F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(s)\le P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}\frac{\mathrm{d}t}{\mathrm{d}s}(s)\end{array}\right.

\left\{\begin{array}{l}\mathrm{0}\le E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)\le E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s)\\ E_{\mathrm{k}\mathrm{i}\mathrm{n}}(S)=E_{\mathrm{k}\mathrm{i}\mathrm{n}}(\mathrm{0})\\ E_{\mathrm{b}}(S)\le E_{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right.（33）

式中：S为单圈赛道长度.目标函数中的\mathrm{d}t/\mathrm{d}s(s)表示时间在空间上的导数，其代表了行驶每段单位距离的时间消耗，通常也被定义为速度v(s)的倒数：

\frac{\mathrm{d}t}{\mathrm{d}s}(s)=\frac{\mathrm{1}}{v(s)}（34）

车辆动能E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)与速度v(s)的关系如下：

E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)=\frac{\mathrm{1}}{\mathrm{2}}mv{(s)}^{\mathrm{2}}（35）

通过式（34）与式（35）即可建立E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)与\mathrm{d}t/\mathrm{d}s(s)的关系，从而计算目标函数的值.

式（33）中的F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(s)为位置域的制动力，其对应于时间域的控制量制动功率P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}(t)；F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)是对应于时间域控制量电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}(t)的等效力；F_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}(E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s))为位置域的总阻力，其可由式（2）推导；F_{\mathrm{p}}(s)为位置域的车辆轮上驱动力，由式（8）进一步推导可得：

\frac{F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)}{v(s)}={\alpha }_{\mathrm{p}}{F}_{\mathrm{p}}{(s)}^{\mathrm{2}}+\frac{F_{\mathrm{p}}(s)}{v(s)}（36）

式（33）中的F_{\mathrm{i}}(s)是对应于时间域电池总功率P_{\mathrm{i}}(t)的等效力，由式（4）进一步推导可得：

\frac{F_{\mathrm{i}}(s)}{v(s)}={\alpha }_{\mathrm{b}}{F}_{\mathrm{b}\mathrm{a}\mathrm{t}}{(s)}^{\mathrm{2}}+\frac{F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)}{v(s)}（37）

前述式（34）~式（37）为位置域各变量间的约束，但这些约束均为非线性等式约束而非凸约束，不可参与凸优化，因此还需对这些约束进行凸化.

对比式（33）与式（14）可发现，式（33）相对于 式（14）所示的时间域优化问题减少了一个状态量，即路程s，故时间域中与路程s相关的协态变量{\lambda }_{\mathrm{3}}不参与后续的凸优化求解，前文1.4.2节已提到{\lambda }_{\mathrm{3}}与最优控制策略函数无直接的关系，因此不计算{\lambda }_{\mathrm{3}}的数值解并不会对后文产生直接影响，{\lambda }_{\mathrm{3}}对{\lambda }_{\mathrm{1}}的影响不在本文的讨论范围内.

2.2 位置域约束凸化

2.2.1 位置域车辆纵向动力学凸约束

前文2.1节中，\mathrm{d}t/\mathrm{d}s(s)与v(s)是凸优化问题中两个不同的优化变量，式（34）是两者间的直接关系约束，通过重新排列与松弛化^［27］，可以得到：

v(s)\frac{\mathrm{d}t}{\mathrm{d}s}(s)\ge \mathrm{1}（38）

进一步推导可得：

{\left(\frac{\mathrm{d}t}{\mathrm{d}s}(s)+v(s)\right)}^{\mathrm{2}}\ge {\mathrm{2}}^{\mathrm{2}}+{\left(\frac{\mathrm{d}t}{\mathrm{d}s}(s)-v(s)\right)}^{\mathrm{2}}（39）

最终即得到一个二阶锥形式的凸约束^［28］：

\frac{\mathrm{d}t}{\mathrm{d}s}(s)+v(s)\ge {‖\begin{array}{c}\mathrm{2}\\ \frac{\mathrm{d}t}{\mathrm{d}s}(s)-v(s)\end{array}‖}_{\mathrm{2}}（40）

在最小化式（33）的目标函数时，最优情况是\mathrm{d}t/\mathrm{d}s(s)尽可能取最小值，因此在问题最优解中一定会满足式（40）的等式情况.当车速v为零时，\mathrm{d}t/\mathrm{d}s会趋近于无穷大，因此式（40）会出现数值问题，但这种情况只发生在车辆起步瞬间或者出现某些事故的情况，本文不考虑这些特殊情况.

前文式（35）可进行与式（38）类似的松弛化，从而得到如下式子：

E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)\ge \frac{\mathrm{1}}{\mathrm{2}}mv{(s)}^{\mathrm{2}}（41）

根据式（38）~式（40）的思想，可得二阶锥形式的凸约束如下所示：

\frac{\mathrm{2}}{m}+E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)\ge {‖\begin{array}{c}\mathrm{2}v(s)\\ \frac{\mathrm{2}}{m}-E_{\mathrm{k}\mathrm{i}\mathrm{n}}(s)\end{array}‖}_{\mathrm{2}}（42）

凸优化问题的目标是最小化单圈时间，因此最优解中的v(s)一定会尽可能取得最大值，从而满足式（42）的等式情况.

2.2.2 位置域车辆动力系统凸约束

对前文式（36）进一步松弛后可得：

\frac{F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)}{v(s)}\ge {\alpha }_{\mathrm{p}}{F}_{\mathrm{p}}{(s)}^{\mathrm{2}}+\frac{F_{\mathrm{p}}(s)}{v(s)}（43）

写为与式（39）类似的形式：

(F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)-F_{\mathrm{p}}(s))\frac{\mathrm{d}t}{\mathrm{d}s}(s)\ge {\alpha }_{\mathrm{p}}{F}_{\mathrm{p}}{(s)}^{\mathrm{2}}（44）

最终即得到一个二阶锥形式的凸约束：

\frac{\mathrm{d}t}{\mathrm{d}s}(s)+F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)-F_{\mathrm{p}}(s)\ge {‖\begin{array}{c}\mathrm{2}\sqrt[]{{\alpha }_{\mathrm{p}}}{F}_{\mathrm{p}}(s)\\ \frac{\mathrm{d}t}{\mathrm{d}s}(s)-F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)+F_{\mathrm{p}}(s)\end{array}‖}_{\mathrm{2}}（45）

式（45）不等式成立的情况对应于部分电能消失的非物理情况，因此在最优解中一定满足式（45）的等式情况.利用前文相同思路对式（37）进行转化可得：

\frac{\mathrm{d}t}{\mathrm{d}s}(s)+F_{\mathrm{i}}(s)-F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)\ge {‖\begin{array}{c}\mathrm{2}\sqrt[]{{\alpha }_{\mathrm{b}}}{F}_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)\\ \frac{\mathrm{d}t}{\mathrm{d}s}(s)-F_{\mathrm{i}}(s)+F_{\mathrm{b}\mathrm{a}\mathrm{t}}(s)\end{array}‖}_{\mathrm{2}}（46）

同样，在最优解中一定满足式（46）的等式情况.

结合2.1~2.2节内容，可得到完整的位置域最小化单圈时间凸优化问题，此问题是完全凸的，可利用现有的二阶锥规划求解器对其进行数值求解，得到协态变量{\lambda }_{\mathrm{1}}与{\lambda }_{\mathrm{2}}的单圈最优变化轨迹数值解，从而得到最优控制量解析解的数值结果.

3 最优控制策略验证

3.1 位置域凸优化问题数值求解方法

前文所述的凸优化问题进行适当离散后便可通过凸优化数值求解器进行求解，此离散优化问题在MATLAB中使用YALMIP^［29］工具箱进行建模，采用MOSEK^［30］求解器进行求解.

在此离散优化问题建模中，路程离散点数量为N=S/\mathrm{\Delta }s，S为单圈赛道长度，路程离散间隔\mathrm{\Delta }s=\mathrm{0.01}\mathrm{m}.两个状态量车辆动能E_{\mathrm{k}\mathrm{i}\mathrm{n}}、电池总能耗E_{\mathrm{b}}，以及两个与动能相关的变量车速v、车速倒数\mathrm{d}t/\mathrm{d}s均是长度为N的一维数组，对两个状态量采用前向欧拉法进行迭代：

x_{k+\mathrm{1}}=x_k+f(u_k)\mathrm{\Delta }s（47）

式中：k=\{\mathrm{1,2},\cdots ,N-\mathrm{2},N-\mathrm{1}\}；f(u_k)为与控制量相关的表达式.两个控制量电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}、制动力F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}，以及其余的变量如车辆轮上驱动力F_{\mathrm{p}}、位置域等效电池总功率F_{\mathrm{i}}、总阻力F_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}等均是长度为N-1的一维数组.在此基础上根据第2节的内容构建数组形式的二阶锥约束与不等式约束.

本文验证所使用的单圈赛道如图5所示.赛道中每个位置s处均有最大限速E_{\mathrm{k}\mathrm{i}\mathrm{n},\mathrm{m}\mathrm{a}\mathrm{x}}(s)，直线区域的限速为一定值且最大，弯道区域的限速设为定值且小于直线区域限速，弯道曲率越大则限速越低.

3.2 关键协态变量数值求解结果

本文在配置有AMD Ryzen 7 6800H处理器和16.0 GB内存的工作站上对某FSAC方程式赛车单圈比赛进行仿真求解.在如图5所示赛道工况下求解得到的协态变量{\lambda }_{\mathrm{1}}、{\lambda }_{\mathrm{2}}单圈最优变化轨迹数值解以及式（29）中的{\lambda }_{\mathrm{1}}临界值{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}、{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}}如图6所示.

由图6可见{\lambda }_{\mathrm{1}}的最优数值解是时变的，{\lambda }_{\mathrm{2}}的最优数值解为常数，此结果与前文1.4.2节的推导结论相匹配.

3.3 最优控制策略数值解验证

图7展示了由高效凸优化求解器计算得到的最优控制量P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}、F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}数值解与{\lambda }_{\mathrm{1}}的关系散点图.前文图4展示的最优控制策略函数解析解P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})、F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})曲线与图7的数据点变化趋势相吻合，最优制动力F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}达到最大值时对应于图4（b）中的最大制动功率P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}.即本文推导的解析形式最优控制策略函数F_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})、F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}({\lambda }_{\mathrm{1}})与凸优化最优解相匹配，其在当前条件下是合理的.

求得关键协态变量{\lambda }_{\mathrm{1}}、{\lambda }_{\mathrm{2}}的数值解后，即可计算P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}的两个临界{\lambda }_{\mathrm{1}}值{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}与{\lambda }_{\mathrm{b}\mathrm{a}\mathrm{t},\mathrm{0}}，在此基础上结合1.4.2节推导出的解析形式最优控制策略函数（27）、式（29）与式（32），即图4，计算赛道每个位置s处最优控制量解析解的数值结果P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}与F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)，F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)对应P_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s).计算流程如表1所示.

图8展示了由表1流程计算得到的最优电池输出功率解析解结果P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)曲线、凸优化求解器数值计算出的P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)曲线以及两者间的误差曲线.由图8（a）可直观地看出两种结果的曲线匹配程度较高.结合图8（a）与图8（b）的误差值与误差率曲线可知，P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)解析解与凸优化结果间的最大误差值为146W，且误差较大值均出现在P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)结果较大的路段，全程大部分路段的误差值小于25 W；最大误差率为6.8%，且仅在少部分路段中出现了4%~6.8%的较大误差率，此部分路段中P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)结果本就较小，因此较小的误差值带来了较大的误差率，全程大部分路段的误差率均小于1%.图8（b）的误差率曲线并非连续，缺失的曲线部分对应的路段中P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)解析解为零，由图8（b）可见这部分路段的误差值均极小.上述提到的误差来自{\lambda }_{\mathrm{1}}、{\lambda }_{\mathrm{2}}的凸优化数值解误差，以及在表1流程中通过式（31）求解P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}(s)时的数值误差，总体误差在可接受范围内.综上，当前工况下的最优电池输出功率解析解结果足够接近最优解.

图9展示了由表1流程计算得到的最优制动力解析解结果F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)曲线、凸优化求解器数值计算出的F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)曲线以及两者间的误差曲线.由图9（a）可直观地看出两种结果的曲线匹配程度较高.结合图9（a）与图9（b）的误差值与误差率曲线可知，F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)解析解与凸优化结果间的最大误差值为660 N，最大误差率很高，但较大的误差值与误差率只出现在全程中的极个别路段，绝大部分路段中误差值与误差率均接近于零.由前文式（27）可知，解析形式的最优制动力F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)曲线类似于阶跃函数，连续两个路段的F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}解析解结果可能会在零至最大值间突变，而凸优化数值解因其数值误差则可能并不会表现出准确的阶跃形式，这便导致了图9（b）中极少数路段的较大误差.图9（b）的误差率曲线同样并非连续，缺失的曲线部分对应的路段中F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}(s)解析解为零，由图可见这部分路段的误差值均极小.由于绝大多数路段的误差值均接近于零，因此总体误差在可接受范围内.综上，当前工况下的最优制动力解析解结果足够接近最优解.

进一步结合图8与图9进行分析，当最优电池输出功率P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}不为0时，最优制动力F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}为0，当P_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{*}}为0时，F_{\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{k}\mathrm{e}}^{\mathrm{*}}不为0且取得最大值.两个控制量的变化趋势都极陡，这是因为本文模型假设控制指令能在最短时间内完成，暂未考虑控制执行器的性能限制.最优控制量的结果反映出车辆总是以尽可能大的加速度进行加减速，这符合最小化单圈比赛时间的要求.综上，本文的控制策略在当前工况下能实现求解器的最优解.

为验证前文最优控制策略解析解的数值结果合理性，本文利用其对式（33）的位置域状态方程进行前向迭代计算，即模拟最优控制策略解析解的实施，最终输出车辆动能E_{\mathrm{k}\mathrm{i}\mathrm{n}}、电池总能耗E_{\mathrm{b}}的变化曲线，并将其与高效凸优化求解器计算的对应结果进行对比分析.

图10展示了由位置域状态方程前向迭代计算出的最优车辆动能E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}}(s)曲线、凸优化求解器数值计算出的E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}}(s)曲线以及两者间的误差曲线.

由图10（a）可直观地看出两种结果的曲线匹配程度较高，且均未超过最大动能限制.由图10（b）的误差值及误差率曲线可知，最大误差值为375 J，结合图10（b）与图9（b）可发现，图10（b）中动能误差曲线每一次大的阶跃均对应于图9（b）中制动力误差的突变，且最大的动能误差对应于最大的制动力误差，由式（33）可知制动力会直接影响车辆动能，即动能误差值曲线的整体变化趋势是符合逻辑的，图8（b）展示的最优电池输出功率误差率在全路段都极低，因此其误差对动能的影响要远小于制动力误差产生的影响.动能误差率在绝大部分路段均低于4%，仅在极少部分路段出现了超过20%的较大误差率，结合图10（a）与图10（b）可知在这部分路段中，动能误差值相对较高但仍是一个较小值，而E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}}(s)本就是一个较小的结果，因此表现出了较大的动能误差率.整体上，由状态方程前向迭代计算得到的最优动能E_{\mathrm{k}\mathrm{i}\mathrm{n}}^{\mathrm{*}}(s)能在较小误差范围内接近凸优化求解器的最优解.

图11展示了由位置域状态方程前向迭代计算出的最优电池总能耗E_{\mathrm{b}}^{\mathrm{*}}(s)曲线、凸优化求解器数值计算出的E_{\mathrm{b}}^{\mathrm{*}}(s)曲线以及两者间的误差曲线.

由图11（a）可直观地看出两种结果的曲线整体匹配程度较高，但在曲线后段出现了较为明显的误差.结合图11（b）的误差值及误差率曲线可知，误差值整体呈增加的趋势，且当部分路段的最优电池总能耗快速增加时，其误差值也迅速增加，最大误差值为6 970 J，由式（33）可知电池总能耗只与电池总功率有关，而电池总功率与电池输出功率有如式（15）所示的直接关系，因此图8（b）所示的最优电池输出功率误差会直接反映在电池总功率中，从而导致最优电池总能耗在前向迭代计算的过程中持续累积误差，最终产生了较大的误差值.虽然最终累积的误差值较大，但由误差率曲线可知路程后段的电池总能耗误差率呈下降趋势，最大误差值对应的误差率仅为1.8%，全程最大误差率为5.9%，路程前段的误差率较大，这是因为路程前段的E_{\mathrm{b}}^{\mathrm{*}}(s)本就较小，较小的误差值导致了较大的误差率.整体上，由状态方程前向迭代计算得到的最优电池总能耗E_{\mathrm{b}}^{\mathrm{*}}(s)能在较小误差范围内接近凸优化求解器的最优解.

进一步结合状态量与控制量进行分析，图10（a）和图11（a）中两个状态量曲线的变化趋势与图8（a）和图9（a）中两个控制量的大小关系相吻合，车辆动能与电池总能耗的快速增加对应于电池输出功率取得局部最大值，车辆动能的快速减小对应于制动力取得最大值，电池总能耗的增加趋势变缓对应于电池输出功率取得较小值，各曲线的数值误差在可接受范围内，并不影响整体的合理性.综合上述可知，本文的控制策略在当前的赛道与车辆条件下具有合理性与最优性.

4 结 论

本文对于特定赛道条件下某无人驾驶纯电动赛车的单圈完赛时间最短问题进行了研究，提出了解析解与数值解相结合的时间最优控制策略.基于车辆的动力学与动力系统模型构建了时间最短的最优控制问题，在时间域内结合非光滑分析与庞特里亚金极小值原理推导出了与关键协态变量相关的解析形式最优控制策略函数；在位置域内将时间最短的最优控制问题建模为二阶锥形式的凸优化问题；通过MATLAB的YALMIP与MOSEK求解工具进行数值求解得到关键协态变量的数值变化轨迹，进一步求得赛车在每个位置的最优控制策略解析解结果，将最优控制策略解析解的实施结果与高效凸优化求解器的计算结果进行对比，验证了该最优控制策略的合理性与有效性.后续需进一步优化电池模型，考虑电池的动态特性对于最优控制策略的影响.此外也应研究控制策略的赛中在线控制可行性，可结合数值解结果与解析解结论分析关键协态变量的实时变化趋势，或结合模型预测控制的相关思想，从而在赛中根据当前状态实时闭环更新车辆在每个位置的最优控制量.

参考文献

[1] 王浩宇， 左志强， 邓宇翔， 等. 无人驾驶车辆运动控制综 述［J］. 控制理论与应用， 2024， 41（9）： 1507-1522．

[2] WANG H Y， ZUO Z Q， DENG Y X，et al. A survey of motion control for unmanned ground vehicles［J］. Control Theory & Applications， 2024， 41（9）： 1507-1522．（in Chinese）

[3] 孙毅文. 分布式驱动车辆复杂赛道时间最优控制研究［D］. 北京： 清华大学，2023．

[4] SUN Y W ．Minimum lap time control for distributed drive vehicle in complex racetrack environment［D］. Beijing：Tsinghua University，2023．（in Chinese）

[5] SCIARRETTA A， DE NUNZIO G， OJEDA L L. Optimal ecodriving control：energy-efficient driving of road vehicles as an optimal control problem［J］. IEEE Control Systems Magazine，2015， 35（5）： 71-90．

[6] DONG H X， ZHUANG W C， CHEN B L， et al ．Predictive energy-efficient driving strategy design of connected electric vehicle among multiple signalized intersections［J］. Transportation Research Part C：Emerging Technologies，2022，137：103595．

[7] DONG H X， ZHUANG W C， CHEN B L，et al ．A comparative study of energy-efficient driving strategy for connected internal combustion engine and electric vehicles at signalized intersections［J］．Applied Energy，2022，310：118524．

[8] KIM Y， FIGUEROA-SANTOS M， PRAKASH N，et al ．Co-optimization of speed trajectory and power management for a fuel-cell/battery electric vehicle［J］．Applied Energy，2020，260：114254．

[9] HERRMANN T， PASSIGATO F， BETZ J，et al ．Minimum race-time planning-strategy for an autonomous electric racecar［C］//2020 IEEE 23rd International Conference on Intelligent Transportation Systems （ITSC）．September 20-23，2020，Rhodes，Greece．IEEE，2020：1-6．

[10] PADILLA G P， WEILAND S， DONKERS M C F ．A global optimal solution to the eco-driving problem［J］．IEEE Control Systems Letters，2018，2（4）：599-604．

[11] HAMEDNIA A， SHARMA N K， MURGOVSKI N，et al ．Computationally efficient algorithm for eco-driving over long look-ahead horizons［J］. IEEE Transactions on Intelligent Transportation Systems，2022，23（7）：6556-6570．

[12] LI Y P， TANG X L， LIN X K，et al ．The role and application of convex modeling and optimization in electrified vehicles［J］．Renewable and Sustainable Energy Reviews，2022，153：111796．

[13] BROERE S， VAN KAMPEN J， SALAZAR M ．Minimum-lap-time control strategies for all-wheel drive electric race cars via convex optimization［C］//2022 European Control Conference （ECC）．July 12-15，2022，London，United Kingdom．IEEE，2022：1204-1211．

[14] BORSBOOM O， FAHDZYANA C A， SALAZAR M ．Time-optimal control strategies for electric race cars with different transmission technologies［C］//2020 IEEE Vehicle Power and Propulsion Conference （VPPC）．November 18 - December 16，2020，Gijon，Spain．IEEE，2021：1-5．

[15] VAN KAMPEN J， HERRMANN T， HOFMAN T，et al ．Optimal endurance race strategies for a fully electric race car under thermal constraints［J］．IEEE Transactions on Control Systems Technology，2024，32（3）：818-833．

[16] DE CASTRO R， TANELLI M， ESTEVES ARAÚJO R，et al ．Minimum-time path-following for highly redundant electric vehicles［J］．IEEE Transactions on Control Systems Technology，2016， 24（2）： 487-501．

[17] SALAZAR M， ELBERT P， EBBESEN S，et al ．Time-optimal control policy for a hybrid electric race car［J］. IEEE Transactions on Control Systems Technology，2017，25（6）：1921-1934．

[18] SALAZAR M， BALERNA C， CHISARI E，et al ．Equivalent lap time minimization strategies for a hybrid electric race car［C］//2018 IEEE Conference on Decision and Control （CDC）. December 17-19， 2018，Miami， USA．IEEE，2019：6125-6131．

[19] SALAZAR M， BALERNA C， ELBERT P， et al ．Real-time control algorithms for a hybrid electric race car using a two-level model predictive control scheme［J］. IEEE Transactions on Vehicular Technology，2017，66（12）： 10911-10922．

[20] CARTIGNIJ C， SALAZAR M ．Time-optimal design and control of electric race cars equipped with multi-speed transmissions［C］//2023 IEEE Vehicle Power and Propulsion Conference （VPPC）. October 24-27，2023，Milan，Italy．IEEE，2024：1-4．

[21] DUHR P， SCHALLER M， ARZILLI L，et al. Time-optimal energy management of the formula 1 power unit with active battery path constraints［C］//2021 European Control Conference （ECC）. June 29-July 2，2021，Delft， Netherlands. IEEE， 2022： 913-920．

[22] BETTIOL P， VINTER R. Dynamic programming［M］//Principles of Dynamic Optimization. Cham：Springer Nature Switzerland，2023：643-767．

[23] LYU S W， YANG L Y， SHEN D L，et al ．Analysis of control behavior in eco-driving speed optimization using Pontryagin’s minimum principle［J］. IEEE Access，2024，12：148893-148903．

[24] DUHR P， CHRISTODOULOU G， BALERNA C，et al ．Time-optimal gearshift and energy management strategies for a hybrid electric race car［J］．Applied Energy，2021，282：115980．

[25] ZHANG P， CHEN H B， MA Z W，et al ．Convex optimization based collision avoidance path planning method for unmanned ground vehicles［C］//2024 36th Chinese Control and Decision Conference （CCDC）．May 25-27，2024，Xi’an，China．IEEE，2024： 844-849．

[26] ZHOU W， CAI X， CHEN Y Q，et al ．Decoding the optimal charge depletion behavior in energy domain for predictive energy management of series plug-in hybrid electric vehicle［J］. Applied Energy，2022，316：119098．

[27] SALAZAR M， DUHR P， BALERNA C，et al ．Minimum lap time control of hybrid electric race cars in qualifying scenarios［J］．IEEE Transactions on Vehicular Technology，2019，68（8）：7296-7308．

[28] WANG Z B ．A survey on convex optimization for guidance and control of vehicular systems［J］．Annual Reviews in Control，2024，57：100957．

[29] KONDA M， HOFMAN T， SALAZAR M ．Energy-optimal design and control of electric powertrains under motor thermal constraints［C］//2022 European Control Conference （ECC）．July 12-15，2022，London，United Kingdom．IEEE，2022：1178-1185．

[30] BORSBOOM O， FAHDZYANA C A， HOFMAN T，et al ．A convex optimization framework for minimum lap time design and control of electric race cars［J］．IEEE Transactions on Vehicular Technology，2021，70（9）：8478-8489．

[31] AKHILESHk K， AMIT D， MANISH T. LMI and YALMIP： Modeling and optimization toolbox in MATLAB［C］//2019 International Conference on VLSI， Communication， and Signal Processing （VCAS）. October 21-23，2019， Prayagraj， India. Springer，2021： 507-515.

[32] MOSEK ApS. The MOSEK optimization toolbox for MATLAB［EB/OL］. （2019-11-27）［2024-1-27］

基金资助

湖南省自然科学基金项目(2022JJ30158)

Natural Science