作为计算机图形学的重要分支，人体运动广泛应用于动画影视、游戏、机器人仿真和生物力学等领域^[1]。基于物理模拟的人体运动，由于能够合成符合自然规律的运动片段，成为近年来计算机动画领域中最活跃的研究方向之一^[2]。在现实中，人所处的环境(摩擦力、支持力、姿态等)实时变化，人具有一种对运动进行控制的天生属性，对环境的变化实时响应，确保运动平稳。模拟上述机制的方法通常包含2类，基于几何的运动生成方法往往需要大量运动数据，通过数据驱动合成；相反，基于物理模拟的方法，通过物理控制器生成运动，不需要大数据支持，且得到的运动片段满足真实环境中的物理规律^[3]。

为了模拟真实的人体物理运动，通常对人体采用抽象建模，将人体转化为多个刚体的链接结构；为了模拟肌肉力矩机制，设计物理控制器，通过对骨骼产生力矩作用，驱动人体物理模型运动^[4-6]；为了保证运动系统稳定，人体物理模型引入了反馈机制，通过对系统状态量的改变，来计算控制器中重要关节的改变量，进而将反馈叠加到原始控制器，确保人体稳定运动，避免摔倒、自碰撞等^[6]。

运动控制器的反馈机制为一种$\mathit{\Gamma :}{\mathit{\boldsymbol{R}}^{\rm{n}}} \to {\mathit{\boldsymbol{R}}^m}$的映射，其中${\mathit{\boldsymbol{R^n}}}$为反馈系统观测状态量的向量差$\delta \mathit{\boldsymbol{s}}$的空间，${\mathit{\boldsymbol{R}}^m}$为控制器改变量$\delta \mathit{\boldsymbol{a}}$的空间，$\mathit{\Gamma }$为一种函数关系映射。每一步仿真需要计算观测状态量的改变，然后反馈并叠加原始控制器，从而保证运动系统稳定。这种反馈控制系统中参数较多，控制复杂，如何方便地求解是一个需要解决的问题。

在对反馈控制策略约束求解过程中，往往涉及多个目标函数，例如根关节惩罚项、姿态惩罚项、能量等，需要对它们进行加权求和以转换为单目标函数求解。各目标函数权重的选取对优化的方向、收敛时间、成功与否至关重要；选择失当则极易进入局部极小值导致无法收敛。实施过程需反复尝试，工作烦重，对用户专业素质要求较高。

物理优化求解中可行解空间很小，子代控制器中大部分不满足约束，从而控制器生成的运动容易导致人体的自碰撞、摔倒。另外，由于理想控制器需要具备鲁棒和稳定特性，且能完成长时间运动，然而这种控制器极难获得，需要设置极大子代个数或长时间优化，所以，直接利用多目标进化算法对物理反馈求解会收敛失败。如何提高落到可行域的个体的命中率且降低迭代次数，是一个难点。

为了解决上述问题，对多目标进化算法改进，利用基于区域密度的多层取优保证种群个体的分布多样性，具备各个方向进化能力；利用种群中失败个体构建“禁选区域”预筛选策略，提高在可行域区间中生成个体的命中率；设计基于“剪枝”的多阶段物理优化流程，舍弃鲁棒性较差的个体，加快算法求解。

早期运动控制大多通过手工设计反馈控制器，利用重心位置、速度或其他关键信息作为观测变量，手动设计反馈项进而控制角色。Yin等人^[7]将根关节的线速度、着地脚与重心的水平距离作为观测项，计算反馈并改变摆动脚的目标姿态，对行走运动进行平衡控制。Lee等人^[8]为了使行走更加稳定，拓展了该反馈，增加对着地脚踝的控制，根据当前序列与参考序列中摆动脚的速度差值设计反馈项，利用逆向运动学改变着地脚的高度，达到平稳着地的效果。手动设计反馈要求人员对运动细节十分熟悉，反馈项权重、不同反馈项间的组合关系需要反复实验，烦杂耗时，对设计人员的专业素质有一定要求。

作为低维物理模型的倒立摆也用于运动控制，生成物理行走。Mordatch等人^[9]首先根据地形计算步态数据，然后利用弹簧倒立摆计算人体重心曲线，进而求解全身运动，以生成适应不同地形、具备转向功能的不同风格运动。赵建军等人^[10]通过二阶倒立摆对人体下肢建模，结合步态信息与地面约束对二阶倒立摆约束求解，同时要求满足运动规律与高层用户需求，进而优化关节力矩，生成更加自然真实的下肢动作。由于低维物理方法仅对下肢物理建模，只能生成与足迹相关的下肢运动，如果不借助外部输入，无法真实还原上肢运动^[2]。

运动数据库的发展面临新的机遇。研究人员尝试借助动作捕捉数据或关键帧指导人体物理运动的生成并进行控制，即前馈控制。Liu等人^[5]尝试利用动作捕捉数据采用基于采样的优化算法进行运动生成，并对走跑、跳跃、翻滚等运动的控制，后期则采用CMA-ES进化算法结合回溯来达到丰富优异的运动效果^[11]。Al Borno等人^[4]则利用关键帧数据，结合滑动窗口技术与进化算法进行优化，可以生成多种不同种类物理运动。此外，Lv等人^[12]首先建立人体运动数据和足底力学信息数据库，利用最大后验概率框架完成数据驱动的逆向动力学求解；然后对输入的运动数据估计足底力和人体内部力矩，完成运动编辑、滤波和运动控制等工作；然而该算法需要提前建立人体运动以及匹配的足底力学信息数据库，不适用于本文问题。前馈控制方法虽然生成多种运动，但是其受限于运动数据，属于开环控制，只能生成与运动数据持续时间相同的运动，后续的运动控制就会失效。

反馈控制来源于自动化控制，随着计算机动画发展，近年来被应用于物理控制。线性反馈作为经典控制策略，由于模型简单易于实现，被用于人体物理运动控制系统。Da Silva等人^[13]将平衡反馈策略转化为二次目标函数，通过二次规划求解获得控制力，生成风格化的物理行走运动。Ding等人^[6]在线性反馈基础上，为了解决反馈矩阵参数众多的问题，将反馈矩阵分解为2个矩阵相乘，对原始反馈矩阵降阶，达到缩减参数个数、降低求解时间的目的。Liu等人^[14]在此反馈控制器基础上，通过合并多个不同种类控制器，实现自动跑酷运动。鲁明等人^[15]改变骨骼质量，在原有反馈控制器基础上采用优化算法生成了面向不同骨架的风格化运动。上述方法大大降低了控制策略设计门槛，可方便地设计反馈控制策略，然而求解过程中目标函数较多，权重设计需要丰富的专业经验，普通用户需要多次试验，任务烦重。

本文采用多刚体与链接关节结构对人体进行物理建模，采用文献[5]的人体骨架模型，身高1.6 m，体重62.5 kg，包含17个刚体、17个关节和50个自由度。人体物理模型参数如图 1、表 1所示。

肌肉对骨骼施加力矩作用，驱动人体运动。为了模拟该机制，本文采用比例微分(PD)控制器，结合目标姿态产生力矩，驱动人体物理模型运动^{[4, 16-18]}，公式为

$ \tau = {k_{\rm{p}}}\left( {{\theta _{\rm{d}}} - \theta } \right) - {k_{\rm{d}}}\dot \theta $

(1)

式中，$\tau $为骨骼所受力矩；$k_\text{p}$、$k_\text{d}$为比例与阻尼系数，与骨骼肌肉属性相关；$\theta $、${\theta _{\rm{d}}}$为当前姿态与目标姿态关节的角度，${\dot \theta }$为当前姿态中关节的角速度。

初始控制器来源于运动数据，选取一定时间的关键帧作为目标姿态，结合时空优化方法求解^{[4-5, 11]}。目标姿态间切换采用有限状态机(FSM)进行转换，通常包含2种切换机制：基于时间机制，间隔一定时间后自动切换；基于事件机制，例如摆动脚接触地面^[19]。走路控制如图 2所示。

初始控制器容易导致角色运动不稳定甚至摔倒，原因是未能针对环境变化作出响应。本文提出一种改进的反馈控制策略，公式为

$ \delta \mathit{\boldsymbol{a}} = {\mathit{\boldsymbol{M}}_{\rm{F}}} \cdot \delta \mathit{\boldsymbol{s}} + b $

(2)

式中，${\mathit{\boldsymbol{M}}_{\rm{F}}}$为反馈矩阵；$\delta \mathit{\boldsymbol{a = a - }}\mathit{\boldsymbol{\overline a}}, \delta \mathit{\boldsymbol{s = s - }}\mathit{\boldsymbol{\overline s}} $，其中$\mathit{\boldsymbol{a}}$、$\mathit{\boldsymbol{\overline a}} $分别为反馈后的控制器和原始控制器，$\mathit{\boldsymbol{s}}$、$\mathit{\boldsymbol{\overline s}} $分别为现实观测量和理想观测量；相比文献[6]的反馈控制策略，本文添加了线性补偿项目$b$，用于抵消原始控制器中的随机误差项。

本文采用计算机物理动画领域中经典的协方差矩阵自适应进化策略(CMA-ES)作为基础算法^{[6, 11, 15, 18, 20]}，结合多目标进化，设计基于剪枝多阶段物理求解框架，对运动反馈问题进行求解。在初始控制器基础上，设计反馈控制策略，首先利用失败、优秀个体信息，采用SVM分类器在解空间中构造禁选区域，对父代产生的候选个体进行预筛选，剔除不合格个体，并进行重采样获得满足禁选区域的个体，随后对所有子代个体进行物理仿真；其次，计算多个目标函数得到适应度值，同时更新SVM分类器；对父代与子代联合集合采用基于密度多层取优算法，获得下一代种群；再次，通过基于剪枝的多阶段物理求解算法判断是否进入下一阶段，循环往复，直至结束；最后将解代入反馈控制系统，生成物理运动序列，算法框架如图 3所示。

多目标优化问题定义为由$D$维决策向量，$N$个目标函数以及$m+n$个约束条件所构成的最优化问题^[21]

$ \begin{array}{l} \min \mathit{\boldsymbol{y}} = f\left( \mathit{\boldsymbol{x}} \right) = \left[ {{f_1}\left( \mathit{\boldsymbol{x}} \right),{f_2}\left( \mathit{\boldsymbol{x}} \right), \cdots ,{f_N}\left( \mathit{\boldsymbol{x}} \right)} \right] \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\\ \;\;\;\;\;\;\;\;\;\;{\rm{s}}.\;{\rm{t}}.\;\;{g_i}\left( \mathit{\boldsymbol{x}} \right) \le 0,\;\;\;\;i = 1,2, \cdots ,m\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{h_j}\left( \mathit{\boldsymbol{x}} \right) = 0,\;\;\;\;j = 1,2, \cdots ,n\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{x}} = {\left[ {{x_1},{x_2}, \cdots ,{x_d}, \cdots ,{x_D}} \right]^{\rm{T}}} \in \mathit{\boldsymbol{X}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_{d\_\min }} \le {x_d} \le {x_{d\_\max }},\;\;\;\;d = 1,2, \cdots ,D \end{array} $

(3)

式中，$\mathit{\boldsymbol{x}}$为$D$维决策变量；$\mathit{\boldsymbol{y}}$、$f(\mathit{\boldsymbol{x}})$为目标函数向量，$f_i(\mathit{\boldsymbol{x}})$为第$i$个目标函数；$\mathit{\boldsymbol{X}}$是$\mathit{\boldsymbol{x}}$组成的决策向量空间，$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$是$\mathit{\boldsymbol{y}}$组成的目标函数空间，$f(\mathit{\boldsymbol{x}})$是决策向量空间$\mathit{\boldsymbol{X}}$到目标函数空间$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$的映射，即$f\left(\mathit{\boldsymbol{x}} \right):\mathit{\boldsymbol{X}} \to \mathit{\boldsymbol{ \boldsymbol{\varOmega}; }}{g_i}\left(\mathit{\boldsymbol{x}} \right)$、${h_j}\left(\mathit{\boldsymbol{x}} \right)$分别为不等式与等式约束。

定义1  Pareto支配关系。对于两个不同的决策向量$\mathit{\boldsymbol{u, v}} \in \mathit{\boldsymbol{X}}$，称$\mathit{\boldsymbol{u}}$支配$\mathit{\boldsymbol{v}}$，记作$\mathit{\boldsymbol{u}} \prec \mathit{\boldsymbol{v}}$，当且仅当$\forall i \in \left[{1, N} \right], {f_i}\left(\mathit{\boldsymbol{u}} \right) \le {f_i}\left(\mathit{\boldsymbol{v}} \right) \wedge \exists j \in \left[{1, N} \right], {f_j}\left(\mathit{\boldsymbol{u}} \right) < {f_j}\left(\mathit{\boldsymbol{v}} \right)$。

定义2  Pareto最优解。一个解$\mathit{\boldsymbol{x^*}}$被定义为Pareto最优解，当且仅当$\neg \exists \mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{X:x}} \prec \mathit{\boldsymbol{x^*}}$。在决策向量空间中，所有Pareto最优解组成的集合称为Pareto最优解集

$ P\left( \mathit{\boldsymbol{X}} \right) = \left\{ {\mathit{\boldsymbol{x}}\left| {\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{X}} \wedge \neg \exists \mathit{\boldsymbol{y}} \in \mathit{\boldsymbol{X}},\mathit{\boldsymbol{y}} \prec \mathit{\boldsymbol{x}}} \right.} \right\} $

(4)

式中，$\mathit{P}\left(\mathit{\boldsymbol{X}} \right)$为求解$\mathit{\boldsymbol{X}} $最优解集函数。Pareto最优解集对应的目标函数组成的集合即为Pareto最优前沿，即$\{ \mathit{f}\left(\mathit{\boldsymbol{x}} \right) = {\left[{{f_1}\left(\mathit{\boldsymbol{x}} \right), {f_2}\left(\mathit{\boldsymbol{x}} \right), \cdots, {f_N}\left(\mathit{\boldsymbol{x}} \right)} \right]^T}\left| {\mathit{\boldsymbol{x}} \in \mathit{P}\left(\mathit{\boldsymbol{X}} \right)} \right|\} $。

定义3  Pareto非支配排序。针对$\mathit{\boldsymbol{X}}$，Pareto最优解集为第0层非支配解集，即${\mathit{\boldsymbol{N}}_0} = P\left(\mathit{\boldsymbol{X}} \right)$；种群$S$去除${\mathit{\boldsymbol{N}}_0}$后的Pareto最优解集记作第1层非支配解集，即${\mathit{\boldsymbol{N}}_1} = P\left({\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{N}}_0}} \right)$；以此类推，第$l$层非支配解集为

$ {\mathit{\boldsymbol{N}}_l} = P\left( {\mathit{\boldsymbol{X}} - \left( {{\mathit{\boldsymbol{N}}_0} \cup \cdots \cup {\mathit{\boldsymbol{N}}_{l - 1}}} \right)} \right),l = 1,2, \cdots $

(5)

式中，${\mathit{\boldsymbol{N}}_{l - 1}}$个体优于${\mathit{\boldsymbol{N}}_{l}}$。

MOEAs的不同层之间相互支配，同层个体之间互不支配。若子代个体数小于同一层个体数目时，需要对其筛选取优，遂引入Hypervolume指标。

解集$\mathit{\boldsymbol{A}}$的Hypervolume指标为$\mathit{\boldsymbol{A}}$中所有个体与参考点$a_\text{ref}$在空间$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$中所构成的超立方体的体积

$ \begin{array}{*{20}{c}} {{S_{{a_{{\rm{ref}}}}}}\left( \mathit{\boldsymbol{A}} \right) = }\\ {\mathit{\Lambda }\left( {\bigcup\limits_{a \in P\left( \mathit{\boldsymbol{A}} \right)} {\left\{ {\left( {{f_1}\left( {a'} \right), \cdots ,{f_N}\left( {a'} \right)} \right)\left| {a \prec a' \prec {a_{{\rm{ref}}}}} \right.} \right\}} } \right)} \end{array} $

(6)

式中，$a_\text{ref}$表示比$\mathit{\boldsymbol{A}}$中所有的点都劣的个体；$\mathit{\Lambda }\left(\cdot \right)$为勒贝格测度^[22]，表示被解集·中的所有点支配但不被$a_\text{ref}$支配所有点组成的超立方体。针对解集$\mathit{\boldsymbol{A}}$中个体$a$的Hypervolume指标(或独立支配区域)为

$ {\Delta _S}\left( {a,\mathit{\boldsymbol{A}}} \right): = {S_{{a_{{\rm{ref}}}}}}\left( \mathit{\boldsymbol{A}} \right) - {S_{{a_{{\rm{ref}}}}}}\left( {\mathit{\boldsymbol{A}}\backslash \left\{ a \right\}} \right) $

(7)

由式(7)可知，个体Hypervolume指标需计算2次整个种群Hypervolume指标。若种群个体数较多，计算量庞大；人体物理运动优化涉及高维目标空间，计算Hypervolume指标是NP-hard问题^[23]，时间复杂度随着目标函数个数增加呈现指数级增长。

对第$g$代种群${\mathit{\boldsymbol{X}}_g}$进行取优，计算下一代的初始种群${\mathit{\boldsymbol{X}}_{g+1}}$时，原则上选择${\mathit{\boldsymbol{X}}_g}$的非支配解集${\mathit{\boldsymbol{N}}_0}$中的个体；若$\left| {{\mathit{\boldsymbol{X}}_{g + 1}}} \right| > \left| {{\mathit{\boldsymbol{N}}_0}} \right|$，则还需要选择${\mathit{\boldsymbol{N}}_1} = P\left({{\mathit{\boldsymbol{X}}_g} - {\mathit{\boldsymbol{N}}_0}} \right)$的优秀个体。若采用Hypervolume指标进行取优，虽然保证第2层个体均匀，但取优过程并未考虑第1层个体分布，容易造成“局部聚集”，难以保证解集的多样性，如图 4所示。

图 4中实心圆表示支配解集${\mathit{\boldsymbol{N}}_0}$，空心圆点表示${\mathit{\boldsymbol{N}}_1}$；由于$\left| {{\mathit{\boldsymbol{X}}_{g + 1}}} \right| > \left| {{\mathit{\boldsymbol{N}}_0}} \right|$，第1层的$a \sim g$都被选择，按照式(7)选择最大独立支配区域，则第2层$h$、$i$被选择；上述计算过程只考虑第2层个体信息，造成了第1层$b$、$c$、$d$、$e$与第2层$h$、$i$之间的聚集效应，第1层中相对稀疏的$a$、$f$、$g$并没有在第2层中被选择。原有方法未能很好地评价每个个体在种群中的分布情况，无法保持解集多样性。

为了解决原有算法多层分布取优时却反多样性且计算量大的问题，本文提出一种基于区域密度的取优算法。首先计算某个个体在一定领域范围内所有个体的距离，然后计算相应的区域密度作为拥挤度算子，最终按照拥挤度算子排名对个体取优。

对于集合$\mathit{\boldsymbol{T}}$中的某一个个体$a$，定义以$a$为中心，半径$r$的超球体作为个体$a$关于$\mathit{\boldsymbol{T}}$的领域集合

$ {\mathit{\boldsymbol{L}}_r}\left( a \right) = \left\{ {j\left| {d\left( {j,a} \right) < r,j \in \mathit{\boldsymbol{T}}} \right.} \right\} $

(8)

式中，$d(j, a)$表示个体$j$、$a$的笛卡儿距离。

个体$a$的区域密度为$a$与${\mathit{\boldsymbol{L}}_r}\left(a \right)$内所有个体之间距离的平均值，即

$ D\left( a \right) = \sum\limits_{i = 1}^{\left| {{\mathit{\boldsymbol{L}}_r}\left( a \right)} \right|} {{l_i}/\left| {{\mathit{\boldsymbol{L}}_r}\left( a \right)} \right|} $

(9)

式中，$\left| {{\mathit{\boldsymbol{L}}_r}\left(a \right)} \right|$为集合${\mathit{\boldsymbol{L}}_r}\left(a \right)$中元素的个数；$l_i$为${\mathit{\boldsymbol{L}}_r}\left(a \right)$中第$i$个元素与$a$之间的距离。

实验过程中发现若将个体区域半径设置过小，容易导致某些区域集合中包含极少数目个体，无法表征周围领域个体分布。反之，若将半径设置过大，则区域密度表示为大范围的分布结果，从而使多个个体结果类似，区分度不大。使用个体之间最大距离的1/2为半径进行计算。

${\mathit{\boldsymbol{X}}_g}^\prime $为第$g$代集合，从中选取$\mu $个个体，构成第$g+1$代种子集合${\mathit{\boldsymbol{X}}_{g + 1}}$。区域密度取优的算法如下：

算法1：基于区域密度的取优算法

参数：$\mu, g$；

输入：${\mathit{\boldsymbol{X}}_g}^\prime $；

输出：${\mathit{\boldsymbol{X}}_{g + 1}}$；

begin

计算${\mathit{\boldsymbol{X}}_g}^\prime $非支配解集${\mathit{\boldsymbol{N}}_0} = P\left({{\mathit{\boldsymbol{X}}_g}^\prime } \right)$

${\mathit{\boldsymbol{N}}_{{\rm{outer}}}} = \emptyset //$候选较优解集合

${\mathit{\boldsymbol{N}}_{{\rm{cur}}}} = {\mathit{\boldsymbol{N}}_0}//$计算区域密度个体的集合

${\rm{while}}\left({\mu > \left({\left| {{\mathit{\boldsymbol{N}}_{{\rm{outer}}}}} \right| + \left| {{\mathit{\boldsymbol{N}}_{{\rm{cur}}}}} \right|} \right)} \right)$

${\mathit{\boldsymbol{N}}_{{\rm{tmp}}}} = P\left({{\mathit{\boldsymbol{X}}_g}^\prime - {\mathit{\boldsymbol{N}}_{{\rm{outer}}}}} \right)$下一层的非支配解集

${\mathit{\boldsymbol{N}}_{{\rm{outer}}}} = {\mathit{\boldsymbol{N}}_{{\rm{outer}}}}\bigcup {{\mathit{\boldsymbol{N}}_{{\rm{cur}}}}} $更新${\mathit{\boldsymbol{N}}_{{\rm{outer}}}}$

${\mathit{\boldsymbol{N}}_{{\rm{cur}}}} = {\mathit{\boldsymbol{N}}_{{\rm{tmp}}}}//$更新${\mathit{\boldsymbol{N}}_{{\rm{cur}}}}$

${\rm{endwhile//while}}\left({\mu > \left({\left| {{\mathit{\boldsymbol{N}}_{{\rm{outer}}}}} \right| + \left| {{\mathit{\boldsymbol{N}}_{{\rm{cur}}}}} \right|} \right)} \right)$

计算${{\mathit{\boldsymbol{N}}_{{\rm{OUTER}}}}}$中元素间最长距离$l_\text{max}$

计算集合${{\mathit{\boldsymbol{N}}_{{\rm{cur}}}}}$中元素$a$的领域集合${\mathit{\boldsymbol{L}}_{{l_{\max }}/2\left(a \right)}}$

对${{\mathit{\boldsymbol{N}}_{{\rm{cur}}}}}$中元素$a$在${\mathit{\boldsymbol{L}}_{{l_{\max }}/2\left(a \right)}}$求取其区域密度$D\left(a \right)$

${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}} \leftarrow $区域密度中排名前$\mu - \left| {{\mathit{\boldsymbol{N}}_{{\rm{outer}}}}} \right|$的元素

${\mathit{\boldsymbol{X}}_{g + 1}} = {\mathit{\boldsymbol{X}}_{{\rm{tmp}}}} \cup {\mathit{\boldsymbol{N}}_{{\rm{outer}}}}$

end

物理运动生成优化问题属于非凸优化且不连续问题^[17]，对不合适的解十分敏感，解稍微出错就会导致“病态”运动结果。为了提高较优个体个数，往往需要人为提高采样变异个体的数目，导致物理求解过程中计算量大，优化速度变慢。

为了降低采样数量，同时提高优异个体的生成率，采用“失败学习”策略进行预筛选，充分利用迭代过程中失败个体信息，构建“失败区域”，以防止之后的个体落到该区域中，进而提高优异个体出现率。相比之前方法，其一物理运动控制是一个复杂过程，存在可行区域与不可行区域，分别对应正常的物理运动和非正常运动，且后者占大部分，优化困难；其二，优化求解为一种迭代过程，迭代中的可行区域始终未发生变化但未知，该方法充分利用之前的较差个体先验信息，通过构建SVM分类器，预测可行区域，指导后续较优个体的生成，提高较优个体生成率，进而加快速度。

一般SVM分类器只能针对线性分类，由于高维空间可行区域属于非线性区域，采用RBF高斯核进行高维空间映射进行分类。高斯核为

$ k\left( {{\mathit{\boldsymbol{x}}_i},{\mathit{\boldsymbol{x}}_j}} \right) = \exp \left( { - \frac{{\left\| {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right\|}}{{2{{\left( {k\sigma } \right)}^2}}}} \right) $

(10)

式中，$k$为松弛因子，$\sigma $为高斯函数标准差。当$k$过小，分类结果将产生过拟合，可能会过滤掉可行解，造成可行区域变小；当$k$过大，该分类器在高次特征上的权重实际上衰减得非常快，所以实际上相当于一个低维的子空间，造成可行区域很大，未起到筛选效果。本文选取经验值0.1。

如果生成的变异个体不满足SVM分类器，则需进行重新采样。如果仅在父代个体周围高斯采样，依然容易超出禁选区域界线。鉴于此，在多次重采样过程中，同时采用韦伯衰减方程临时缩小优化步长，保证新个体远离禁选区域，公式为

$ {f_{{\rm{weibull}}}}\left( t \right) = {{\rm{e}}^{ - {{\left( {t/L} \right)}^k} \times \ln \left( 2 \right)}} $

(11)

式中，$L$表示半衰周期，选取10；$k$为衰减的形状参数，选取较快的衰减，设置为2；$t$为尝试次数，避免步长太小，最多容许15次尝试。

为了保证控制器的鲁棒性，本文要求控制器至少可以完成30个周期的完整运动。由于人体物理求解属于时序相关问题，若直接采用MO-CMA及其改进算法求解，会因为低周期个体数过多，容易导致收敛失败，如图 5所示。

图 5为2个目标函数，20个子代个体时多目标优化示意图，其中$t$为控制器的个体可以完成物理运动的周期数，$f_1$为时间惩罚函数，$f_2$为自定义惩罚函数，比如姿态惩罚项。可以看出较高周期的控制器比较少。一方面由于较高周期控制器需要由较低周期控制器进化而来，变异后的优秀解较少；另一方面较低周期控制器会因为其他目标函数(姿态、能量等)值较小的原因，从而形成非支配解，依然在子代中占据大部分。倘若较低周期控制器处理不当，会导致整个种群丧失朝稳定、鲁棒方向进化的能力，直至收敛失败。

针对MO-CMA算法在物理反馈求解中无法收敛的问题，提出一种基于剪枝思想的多阶段物理求解流程，核心思想是在优化一段时间后，且较高周期个体占据一定份额，对较低周期个体进行舍弃，从而保证整个优化流程不在较低个体之间徘徊。如图 6所示。

图 6中，圆圈中的数字表示个体控制器可以完成的周期数，优化窗口表示该阶段正在该窗口周期内进行优化。倘若进化过程中出现超出该窗口的个体，则将其认为优秀个体，反之认为较差个体。本文认为一旦优秀个体占总体的20%后，则表示该种群可以进入下一阶段：当前工作窗口往后移动2个周期，并舍弃掉不满足新窗口的较差个体。该方法一方面保证了较高控制器个体占比，使种群总体朝长时间进化；另一方面由于较低个体占主体，也保证种群具备朝其他目标优化的能力。

为了保证个体产生更多满足并且超越当前窗口的个体，本文采用了基于SVM的预筛选策略保证种群中生成优秀个体的能力。现有初始种群${\mathit{\boldsymbol{X}}_{{\rm{init}}}}$，物理仿真周期数$l_\text{loop}$，工作窗口为$\left[{{w_{{\rm{lower}}}}, {w_{{\rm{upper}}}}} \right]$，优秀个体占比为$r_\text{better}$，物理反馈优化求解流程如下：

算法2：物理反馈优化求解流程

参数：${l_{{\rm{loop}}}}, \left[{{w_{{\rm{lower}}}}, {w_{{\rm{upper}}}}} \right], {r_{{\rm{better}}}}$；

输入：${\mathit{\boldsymbol{X}}_{{\rm{init}}}}$；

输出：${\mathit{\boldsymbol{X}}_{{\rm{out}}}}$；

begin

$g = 0, {w_{{\rm{lower}}}} = 0, {w_{{\rm{upper}}}} = 3$

${\mathit{\boldsymbol{X}}_g} = {\mathit{\boldsymbol{X}}_{{\rm{inin}}}}$第$g$代种群

${\rm{while}}\left({{w_{{\rm{lower}}}} < {l_{{\rm{loop}}}}} \right)$

${\rm{while}}\left({{r_{{\rm{better}}}} < 0.2} \right)$

在${\mathit{\boldsymbol{X}}_g}$中采用SVM预筛选产生后代${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}}$

对后代${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}}$进行物理仿真

对${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}}$筛选出较差个体集合${\mathit{\boldsymbol{X}}_{{\rm{neg}}}}$

${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}} = {\mathit{\boldsymbol{X}}_{{\rm{tmp}}}} - {\mathit{\boldsymbol{X}}_{{\rm{neg}}}}$

利用${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}}$, ${\mathit{\boldsymbol{X}}_{{\rm{neg}}}}$更新SVM分类器

${\mathit{\boldsymbol{X}}_g} \leftarrow $算法1对集合${\mathit{\boldsymbol{X}}_{{\rm{tmp}}}} \cup {\mathit{\boldsymbol{X}}_\mathit{g}}$取优

$g=g+1$

在${\mathit{\boldsymbol{X}}_g}$中计算较优个体比率${{r_{{\rm{better}}}}}$

${\rm{endwhile//while}}\left({{r_{{\rm{better}}}} < 0.2} \right)$

${w_{{\rm{lower}}}} = {w_{{\rm{lower}}}} + 2, {w_{{\rm{upper}}}} = {w_{{\rm{upper}}}} + 2$

${\rm{endwhile//while}}\left({{w_{{\rm{lower}}}} < {l_{{\rm{loop}}}}} \right)$

本文对人体进行物理建模，将肢体抽象为刚体，关节抽象为铰链关节。对物理人模型进行仿真过程，需要满足机器人动力学方程

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{M}}\left( \mathit{\Theta } \right)\mathit{\ddot \Theta } + \mathit{\boldsymbol{V}}\left( {\mathit{\dot \Theta },\mathit{\Theta }} \right) + }\\ {{\mathit{\boldsymbol{J}}^{\rm{T}}}\left( \mathit{\Theta } \right) \times \mathit{\boldsymbol{F}} + \mathit{\boldsymbol{G}}\left( \mathit{\Theta } \right)} \end{array} $

(12)

式中，$\mathit{\boldsymbol{M}}$、$\mathit{\boldsymbol{V}}$、$\mathit{\boldsymbol{G}}$、$\mathit{\boldsymbol{\tau }}$分别表示动力学系统中惯量矩阵、非惯性项、重力项和空间广义力；$\mathit{\boldsymbol{F}}$为外界作用力，如摩擦力、支持力等；$\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$表示广义坐标向量，${\mathit{\boldsymbol{ \boldsymbol{\dot \varTheta} }}}$和${\mathit{\boldsymbol{ \boldsymbol{\ddot \varTheta} }}}$分别表示其一阶导数和二阶导数；$\mathit{\boldsymbol{J}}$为广义坐标系到世界坐标系下的雅可比矩阵。

本文采用全物理仿真，实时判断物理模型与地面的摩擦力和支持力，采用库仑摩擦锥模型^[24]

$ \mathit{\boldsymbol{F}} = \mathit{\boldsymbol{\mu V}} $

(13)

式中，$\mathit{\boldsymbol{F}}$为地面作用于物理模型的力，$\mathit{\boldsymbol{V}}$为4组摩擦锥线性基，$\mathit{\boldsymbol{\mu }}$是4组摩擦锥的加权系数组成的向量。

在进行迭代训练控制器过程中，需要人为设置一些禁选区域约束条件。当子代生成的运动不满足该类约束，则表示该个体处于禁选区域。

在优化窗口范围内对个体进行求解，期望个体能完成窗口规定的最小周期的运动，即

$ \begin{array}{*{20}{c}} {{c_{{\rm{time}}}}\left( i \right) = }\\ {\left\{ {\begin{array}{*{20}{c}} 1&{\left( {{t_e} - {w_{{\rm{lower}}}}\left( i \right)} \right) \le 0}\\ 0&{\left( {{t_e} - {w_{{\rm{lower}}}}\left( i \right)} \right) > 0} \end{array}} \right.,\;\;\;i = 1,2,3,4} \end{array} $

(14)

式中，${c_{{\rm{time}}}}\left(i \right)$表示第$i$个优化阶段时间约束，$t_e$表示角色在该阶段可以完成的周期数，${w_{{\rm{lower}}}}\left(i \right)$表示第$i$个阶段的优化窗口最低需要完成的周期数。满足该约束表示角色能正常完成该优化过程；否则表示中途因跌倒或者弹飞等原因中止优化过程。

物理模型由多个刚体通过关节链接而成，其运动过程中不允许发生自碰撞，故在反馈学习迭代过程中设计自相交约束

$ {c_{{\rm{selfcolide}}}} = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 1\\ 0 \end{array}&\begin{array}{l} 未发生自碰撞\\ 发生自碰撞 \end{array} \end{array}} \right. $

(15)

为了保证控制器的鲁棒性，在训练个体时增加一定的外部干扰力。然而，外部力施加到物理角色后，会导致人体姿态变形，与参考序列相差变大，影响整个区间的目标函数值。这种应激反应区间不属于目标函数的考查区间。因此，不考虑施加力后面的3个周期的目标函数值，允许角色在该段时间自动恢复，如图 7所示。

为了衡量物理模拟运动的效果，设计以下目标函数:

1) 根关节惩罚项目。对于普通运动，本文希望物理仿真运动与参考序列的根关节朝向尽可能相似。设计根关节惩罚项

$ {E_{{\rm{root}}}} = \frac{1}{T}\int {{D_q}\left( {{q_{{\rm{root}}}}\left( t \right),{q_{{\rm{root\_ref}}}}\left( t \right)} \right){\rm{d}}t} $

(16)

式中，${D_q}\left({{q_1}, {q_2}} \right)$表示四元数$q_1$、$q_2$之间的距离，${q_{{\rm{root}}}}\left(t \right)$、${q_{{\rm{root\_ref}}}}\left(t \right)$表示$t$时刻物理仿真与参考序列的根关节朝向四元数。

2) 姿态惩罚项。为了提高物理仿真运动效果，保证物理仿真运动各个关节的姿态与参考序列相似。设计姿态惩罚项

$ {E_{{\rm{pose}}}} = \frac{1}{{T \times {j_{{\rm{num}}}}}}\int {\sum\limits_{j = 1}^{{j_{{\rm{num}}}}} {{w_j}\left( {{D_q}\left( {{q_j}\left( t \right),{q_{j{\rm{\_ref}}}}\left( t \right)} \right)} \right){\rm{d}}t} } $

(17)

式中，${q_j}\left(t \right)$、${q_{j\_{\rm{ref}}}}\left(t \right)$分别表示物理仿真与参考序列中第$j$个关节四元数朝向，$w_j$为每个关节在姿态惩罚中的权重，本文认为所有关节是等价的，关节权重为1。

3) 运动对称惩罚项。对于行走类运动，其在左脚或者右脚着地时刻，运动呈现镜面对称特性。因此，为了使行走类运动看起来左右对称均衡，需要进行对称项目惩罚，即

$ \begin{array}{*{20}{c}} {{E_{{\rm{sym}}}} = \frac{1}{{2 \times p}}\sum\limits_{i = 1}^{2 \times p} {\left( {\frac{1}{{{j_{{\rm{num}}}}}}\sum\limits_{j = 1}^{{j_{{\rm{num}}}}} {{D_q}\left( {\mathit{\boldsymbol{Q}}\left( {i,j} \right),\overline {\mathit{\boldsymbol{Q}}\left( {i - 1,j} \right)} } \right)} + } \right.} }\\ {\left. {{D_q}\left( {{q_{{\rm{root}}}}\left( i \right),\overline {{q_{{\rm{root}}}}\left( {i - 1} \right)} } \right)} \right)} \end{array} $

(18)

式中，$p$为周期个数，$j_\text{num}$表示关节的总个数，$\mathit{\boldsymbol{Q}}\left({i, j} \right)$、$\mathit{\boldsymbol{\overline Q}} \left({i, j} \right)$、${q_{{\rm{root}}}}\left(i \right)$、${\overline q _{{\rm{root}}}}\left(i \right)$分别表示第$i$次脚着地时物理模型第$j$个关节的朝向、镜面朝向以及根关节朝向、根关节镜面朝向。

4) 能量惩罚项目。人体在进化过程中，为了存储更多体力保护自己，潜意识地遵循着能量消耗最小或接近最小的准则进行运动。依据上述准则，设计能量惩罚项

$ {E_{{\rm{energy}}}} = \frac{1}{T}\int {\sum\limits_{r = 1}^{{r_{{\rm{num}}}}} {\left( {\frac{1}{2}{\mathit{\boldsymbol{I}}_r}\mathit{\boldsymbol{w}}_r^2 + {m_r}g{h_r}} \right){\rm{d}}t} } $

(19)

式中，${r_{{\rm{num}}}}$为物理模型中刚体数目，${\mathit{\boldsymbol{I}}_r}$为第$r$个刚体的惯量，${\mathit{\boldsymbol{w}}_r}$为第$r$个刚体绕父关节转动的角速度，$m_r$、$h_r$分别为第$r$个刚体质心质量与高度。本项目考查在一段时间内所有刚体的转动机械能的大小。

实验硬件环境为Intel Xeon CPU E5-2609 v3处理器，8 GB内存，GeForce Titan X显卡。软件中物理底层采用Open Dynamic Engine (ver 0.14)，物理仿真步长为0.002 5 s。

本文解决面向反馈控制器的优化问题，分别提出了改进的多目标进化算法与面向物理反馈控制优化求解流程。分别针对这两部分进行实验与分析。

为了衡量本文提出的改进多目标进化算法的优劣，采用业内具有代表性的测试函数^[25-27]，并设计相应的约束区域，分别采用MO-CMA与本文算法进行实验对比。测试函数和相关参数如表 2所示。

为了衡量多目标进化算法解集的收敛与分布多样性好坏，本文引入反转世代距离(IGD)^[28]指标，用于表征多目标进化算法的解集与真实的Pareto前沿之间的距离。其值越小，表示与真实的Pareto前沿越相似，分布越均匀；反之，则越远，分布越分散。公式为

$ {D_{{\rm{IGD}}}} = \frac{1}{{\left| {{\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}} \right|}}\sum\limits_{i = 1}^{\left| {{\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}} \right|} {{D_{{\rm{mindist}}}}\left( i \right)} $

(20)

式中，${\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}$为真实Pareto解，${D_{{\rm{mindist}}}}\left(i \right) = \mathop {\min }\limits_{j = 1}^{\left| {{\mathit{\boldsymbol{X}}_{{\rm{cur}}}}} \right|} \sqrt {\sum\limits_{m = 1}^N {{{\left({\frac{{{f_m}\left({\mathit{\boldsymbol{x}}_i^{_{{\rm{pareto}}}}} \right) - {f_m}\left({\mathit{\boldsymbol{x}}_j^{{\rm{cur}}}} \right)}}{{f_m^{\max } - f_m^{\min }}}} \right)}^2}} } $表示在${\mathit{\boldsymbol{X}}_{{\rm{cur}}}}$中与${\mathit{\boldsymbol{x}}_i^{_{{\rm{pareto}}}}}$的最小归一化距离，${{\mathit{\boldsymbol{X}}_{{\rm{cur}}}}}$为当前解集，$N$为目标函数的个数，${\mathit{\boldsymbol{x}}_i^{_{{\rm{pareto}}}}}$、${\mathit{\boldsymbol{x}}_j^{{\rm{cur}}}}$分别为${\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}$的第$i$个个体与${\mathit{\boldsymbol{X}}_{{\rm{cur}}}}$的第$j$个个体。${f_m^{\max }}$、${f_m^{\min }}$分别为第$m$个目标函数值的最大、最小值。由于带约束区域目标函数的${\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}$并不可知，本文默认2 000次迭代后的解集近似为${\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}$。

如何更多地生成满足约束区域子代个体也是本部分算法考查指标之一，命中率指标$R_\text{hit}$为

$ {R_{{\rm{hit}}}} = \frac{{{N_{{\rm{statify}}}}}}{{{N_{{\rm{total}}}}}} $

(21)

式中，$N_\text{statify}$为满足约束区域的子代个体数，$N_\text{total}$为所有子代个体数。分别采用MO-CMA与本文算法对3种测试函数进行实验，计算IGD及命中率，如图 8所示。

由图 8可以看出，在迭代次数相同的情况下，本文方法的IGD指标比原始方法小，表明本文方法的解集更加接近于${\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}$，且分布更加均匀；在迭代次数相同的情况下，本文算法的命中率比原始方法高，说明本文算法可以生成更多的满足约束区域的子代，从而更容易产生优秀子代。另外，2种方法的IGD指标总体上随着进化代数的增加而变小，表明方法总体是朝向${\mathit{\boldsymbol{X}}_{{\rm{pareto}}}}$方向进化；算法的命中率刚开始很低，说明初始种群并未在约束区域中，处于寻找约束区域阶段，待寻找到可行区域后，命中率就保持一定数值。

然而，图 8(c)中DTLZ4函数的原始方法的IGD指标在100次后突然上升，并且图 8(f)命中率也随之上升，通过实验分析发现是由进入局部解后分布变得聚集导致的。如图 9所示。

由图 9可以看出，2种算法在80次时分布相差不大，MO-CMA的解集在左半部分分布较少，大部分集中在右边区域。然而随着进化代数增加，在200次时，MO-CMA算法由于进入局部解，分布逐渐偏向于一条直线，分布性变差，从而导致IGD指标突然变大；对于命中率来说，迭代后期，种群的大部分个体处于局部极小值中，生成的新的个体也被局限在该区域，加之步长随着代数增加而减少，从而导致命中率变高。另一方面，本文算法一直保持良好的分布性，解集比之前均匀。

通过实验发现，直接采用MO-CMA算法并不能解决人体物理优化求解问题，参见4.1节，较低周期的个体占据子代个体中的大部分，使种群丧失朝长周期进化的能力，进而容易导致收敛失败。因而，物理反馈的优化求解在改进的MO-CMA算法基础上，采用基于剪枝多阶段物理优化流程进行求解。本部分针对走路、跑步和翻滚运动，不仅涉及经典的行走控制器，镜面对称运动，并且也对非对称运动控制器进行求解。物理运动及实验相关参数如表 3所示。

走路与跑步运动由于具有镜面对称特性，其控制器可以利用该性质将优化变量减半。在实验分析中发现目标函数的个数直接影响到求解空间的大小。若将目标函数的个数设置过多，子代个数必须随之成倍增加，才能保证优化窗口顺利后移，势必会造成求解难度增加。因此，设置3个目标函数，其中包含：1)周期相关目标函数，保证朝着稳定和鲁棒性进化；2)姿态相关目标函数，保证与参考序列尽可能相似，若该运动涉及镜面对称，加入对称惩罚项；3)能量相关目标函数，保证运动尽可能“省力”。

在约束条件中，所有运动需满足不自碰撞约束；约束求解是按照优化窗口进行的，所有当前个体控制下的物理运动，完成运动周期数目不能小于当前窗口的最小数目，否则个体予以丢弃；对于行走运动，则满足稳定约束，不能摔倒，然而此限制条件并不适用于翻滚，为了保证翻滚正常求解，引入根关节约束，需要保证根关节与参考序列的根关节在一定范围内，保证角色整体朝向与参考序列基本一致，即

$ {D_q}\left( {{q_{{\rm{root}}}},{q_{{\rm{root\_ref}}}}} \right) < 1 $

(22)

利用算法生成物理反馈控制的运动效果如图 10所示，其中红色的直线为物理模型与地面的受力，具体动画效果请参考论文附带视频。

控制器的求解涉及长周期的学习训练过程，单纯利用MO-CMA算法无法对其进行求解，因此将MO-CMA算法与基于剪枝的多阶优化相结合，并与本文算法进行对比。其进化过程中产生个体的命中率及目标函数值分别如图 11和图 12所示。

由图 11可以看出，本文方法命中率高于另一种方法，保证每一轮迭代产生更多优秀子代，进而更快地达到当前窗口往后移动的条件，进入下一阶段优化，加快完成物理优化求解。从图 11看出，本文算法的迭代次数更少，时间性能更优。单看某一种算法，其命中率的总体趋势是越来越小，原因是随着优化窗口往后移动，对个体的要求也随之增加，个体控制器需具有完成更长周期的能力，满足此条件的子代个体越来越少，命中率也随之降低。

本文分别对2种算法中所有子代个体目标函数值取平均，计算结果如图 12所示。当迭代结束后，所有个体都能完成规定周期的运动，其第1个目标函数值是一致的，遂不进行对比。从图 12中可以看出，本文算法的最终子代个体目标函数值低于其他算法。由于第3个子目标函数为能量函数，其目标值比较大，为了方便画图，遂对其进行放缩。

为了让控制器更加稳定、鲁棒，在优化训练过程中，添加了部分外部力干扰。人体物理模型在受力后控制器生成的反馈及自我恢复运动序列如图 13所示。

结合算法描述和实验结果与分析，现将本文提出的面向反馈控制器的多目标求解方法总结如下：

针对多目标算法中选取优秀个体可能分布不均匀、优秀个体命中率低的问题，采取区域密度多层取优、禁选区域预筛选求解，达到均匀选取子代的目的，同时提高了优秀个体命中率。

针对人体物理运动求解中多个目标函数权重选取复杂困难的问题，采用多目标进化算法思想，将多个目标等价同时进行约束求解，省去大量实验调参过程；同时提出了基于剪枝的多阶段物理优化流程，通过优化窗口舍弃较差个体，加快物理运动优化求解速度。

然而，本文算法也存在一定的局限性。本文需人工选取控制器的观测项和反馈作用关节，对普通用户门槛高。另外，求取的控制器一次完成一种运动，并不能完成由多种控制器控制生成的长时间运动，不具备控制器切换功能。如何自动选取反馈观测项和作用关节，以及控制器之间的切换将是接下来的工作。