第5章
CHAPTER 5


机械手迭代学习控制





5.1迭代学习控制的数学基础
迭代学习控制是通过迭代修正达到某种控制目标的改善，它的算法较为简单，且能在给定的时间范围内实现未知对象实际运行轨迹以高精度跟踪给定期望轨迹，且不依赖系统的精确数学模型，因而一经推出，就在控制领域得到了广泛的运用。
下面介绍学习控制的稳定性和收敛性分析时用到的基本数学知识［1，2］。
5.1.1矩阵的迹及初等性质
定义设A是n阶方阵，则称A的主对角元素的和为A的迹，记作tr(A)。即若

A=a11a12…a1n
a21a22…a2n
︙︙︙
ai1ai2aiiain
︙︙︙
an1an2…ann（5.1）

则

tr(A)=∑ni=1aii（5.2）

设A、B都是n阶方阵，λ、μ为任意复数，则矩阵的迹具有如下性质： 
（1） tr(λA+μB)=λtr(A)+μtr(B)。
（2） tr(A)=tr(AT)。
（3） 若A∈Cm×n，B∈Cn×m，则tr(AB)=tr(BA)。
5.1.2向量范数和矩阵范数
1. 向量范数

任取x∈Cn，且x=(ξ1ξ2…ξn)T，可定义

‖x‖1=∑ni=1ξi （5.3）

‖x‖2=∑ni=1ξi21/2（5.4）

‖x‖∞=max1≤i≤nξi（5.5）

上述三个范数都是Cn中的向量范数，分别称为1范数，2范数和∞范数，这三个范数实际上都是p范数的特殊情形。
p范数定义如下： 

‖x‖p=∑ni=1ξip1/p，1≤p<+∞（5.6）






2. 矩阵范数
定义： 若对任意矩阵A∈Cm×n，都有实数‖A‖与之对应，且满足下面的范数公理： 
(1) 正定性： ‖A‖≥0，当且仅当A=0时‖A‖=0。
(2) 齐次性： 对任何λ∈C，‖λA‖=λ‖A‖。
(3) 三角不等式： 对任何A,B∈Cm×n，有
‖A+B‖≤‖A‖+‖B‖

则称这个实数‖A‖为矩阵A的范数。

‖A‖V1=def∑mj=1∑ni=1aij（5.7）

‖A‖V∞=defmaxi,jaij（契比雪夫范数）（5.8）

‖A‖Vp=def∑mj=1∑ni=1aijp1/p,1≤p≤+∞（5.9）

当p=2时，称‖A‖V2=‖A‖F=def∑mj=1∑ni=1aij21/2为A的Frobenius范数，简称F范数，是最常用的范数之一。它就是酉矩阵Cm×n中的内积A|B=tr(BHA)所诱导的范数： 

‖A‖2F=(A|A)=tr(
BHA)=∑mj=1∑ni=1aij2（5.10）

F范数具有下列良好的性质： 
性质1设A∈Cm×n，对酉矩阵U∈Cm×m，V∈Cn×n，恒有

‖A‖F=‖UA‖F=‖AV‖F=‖UAV‖F（5.11）

性质2设A∈Cm×n，B∈Cn×l，则有

‖A‖F≤‖A‖F‖B‖F（5.12）

性质3在矩阵空间Cn×n上的任意实函数，记为‖·‖，如果对所有的A,B∈Cn×n,λ∈C，都满足以下4个条件： 
(1) ‖A‖≥0，当且仅当A=0时，有‖A‖=0。
(2) ‖λA‖=λ‖A‖。
(3) ‖A+B‖≤‖A‖+‖B‖。
(4) ‖AB‖≤‖A‖‖B‖。

则称‖·‖为相容的矩阵范数，或简称矩阵范数。显然，矩阵的F范数是一种相容的矩阵范数。
5.2迭代学习控制方法介绍
迭代学习控制（iterative learning control,ILC）是智能控制中具有严格数学描述的一个分支。1984 年，Arimoto［1］等提出了迭代学习控制的概念，该控制方法适合于具有重复运动性质的被控对象，它不依赖于系统的精确数学模型，能以非常简单的方式处理不确定度相当高的非线性强耦合动态系统。目前,迭代学习控制在学习算法、收敛性、鲁棒性、学习速度及工程应用研究上取得了巨大的进展。
近年来，迭代学习控制理论和应用在国外得到快速发展，取得了许多成果。在国内，迭代学习控制理论也得到广泛的重视，有许多重要著作出版［2~5］，发表了许多综述性论文［6~9］。
5.2.1迭代学习控制基本原理
设被控对象的动态过程为

x·(t)=f(x(t),u(t),t)，y(t)=g(x(t),u(t),t)(5.13)

其中，x∈Rn、y∈Rm、u∈Rr分别为系统的状态、输出和输入变量，f(·)、g(·) 为适当维数的向量函数，其结构与参数均未知。若期望控制ud(t)存在,则迭代学习控制的目标为： 给定期望输出yd(t) 和每次运行的初始状态xk(0),要求在给定的时间t∈［0,T］内,按照一定的学习控制算法通过多次重复的运行,使控制输入uk(t)→ud(t)，而系统输出yk(t)→yd(t)。第k次运行时，式(5.13) 表示为

x·k(t)=f(xk(t),uk(t),t),yk(t)=g(xk(t),uk(t),t)(5.14)

跟踪误差为

ek(t)=yd(t)-yk(t)(5.15)

迭代学习控制可分为以下开环学习和闭环学习两种方法： 
(1) 开环学习控制的方法是： 第k+1次的控制等于第k次控制再加上第k次输出误差的校正项，即

uk+1(t)=L(uk(t),ek(t))(5.16)

(2) 闭环学习控制的方法是： 取第k+1次运行的误差作为学习的修正项，即

uk+1(t)=L(uk(t),ek+1(t))(5.17)

其中，L为线性或非线性算子。
5.2.2基本的迭代学习控制算法
Arimoto等首先给出了线性时变连续系统的D型迭代学习控制律［1］

uk+1(t)=uk(t)+Γe·k(t)(5.18)

其中，Γ为常数增益矩阵。在D型算法的基础上，相继出现了P型、PI型、PD型迭代学习控制律。从一般意义来看它们都是PID型迭代学习控制律的特殊形式，PID迭代学习控制律表示为

uk+1(t)=uk(t)+Γe·k(t)+Φek(t)+Ψ∫t0ek(τ)dτ(5.19)

其中，Γ、Φ、Ψ为学习增益矩阵。算法中的误差信息使用ek(t)称为开环迭代学习控制，如果使用ek+1(t) 则称为闭环迭代学习控制，如果同时使用ek(t)和ek+1(t)则称为开闭环迭代学习控制。
此外，还有高阶迭代学习控制算法、最优迭代学习控制算法、遗忘因子迭代学习控制算法和反馈前馈迭代学习控制算法等。
5.2.3迭代学习控制主要分析方法
学习算法的收敛性分析是迭代学习控制的核心问题，这方面的研究成果很丰富。
1. 基本的收敛性分析方法
对于如下线性离散系统： 

x(t+1)=Ax(t)+Bu(t)

y(t)=Cx(t)(5.20)

迭代学习控制算法为

uk+1(t)=uk(t)+Γek(t+1)(5.21)

针对学习算法式(5.21) 的收敛性，有以下两种分析方法： 
（1) 压缩映射方法： 即系统要求满足全局Lipschitz条件和相同的初始条件,如果‖I-CBΓ‖<1，则有

‖ek+1‖=‖(I-CBΓ)ek‖<‖I-CBΓ‖‖ek‖<‖ek‖(5.22)

此时算法是单调收敛的。该方法依赖于范数的选择，常用的有l1范数、l2范数、l∞范数及λ范数。在收敛性证明过程中常用到BellmanGronwall引理。
（2) 谱半径条件法： 如果谱半径ρ满足ρ(I-CBΓ)≤ρ<1,则有

limk→∞‖ek‖= limk→∞‖(I-CBΓ)ek-1‖= limk→∞ρ(I-CBΓ)k‖e0‖(5.23)

即limk→∞‖ek‖=0。
2. 基于2D理论的分析方法
迭代学习控制系统的学习是按两个相互独立的方向进行： 时间轴方向和迭代次数轴方向，因此迭代学习过程本质上是二维系统，可利用成熟的2D系统理论系统地研究和分析时间域的稳定性和迭代次数域的收敛性问题。2D系统的稳定性理论为迭代学习控制的收敛性证明提供了一种非常有效的方法，2D 系统理论中的Roesser模型成为迭代学习控制中最基本的分析模型。
3. 基于Lyapunov直接法的设计方法
Lyapunov直接法已广泛用于非线性动态系统的控制器设计和分析中，在研究非线性不确定系统时，该方法是最重要的应用工具之一。受Lyapunov直接法的启发，在时间域和迭代域能量函数的概念得到研究，它为学习控制在迭代域设计和收敛性分析方面提供了一种新的研究方法。
在迭代域能量函数的迭代学习控制方法基础上，发展了鲁棒和自适应迭代学习控制，可解决具有参数或非参数不确定性非线性系统控制器的设计。近年来反映时间域和迭代域系统能量的组合能量函数方法也应用于迭代学习控制，它可保证在迭代域跟踪误差的渐进收敛以及在时间域具有有界和逐点跟踪的动态特性,并且控制输入在整个迭代区间内是范数收敛的，适用于一类不具有全局Lipschitz 条件的非线性系统。通过能量函数的方法，许多新的控制方法，如反演设计和非线性优化方法都作为系统设计工具应用到迭代学习控制中。
此外，还有最优化分析方法、频域分析法等分析方法。
5.2.4迭代学习控制的关键技术
1. 学习算法的稳定性和收敛性

稳定性与收敛性是研究当学习律与被控系统满足什么条件时，迭代学习控制过程才是稳定收敛的。算法的稳定性保证了随着学习次数的增加，控制系统不发散，但是对于学习控制系统而言，仅仅稳定是没有实际意义的，只有使学习过程收敛到真值，才能保证得到的控制为某种意义下最优的控制。收敛是对学习控制的最基本要求，多数学者在提出新的学习律的同时，基于被控对象的一些假设，给出了收敛的条件。例如，Arimoto在最初提出PID型学习控制律时，仅针对线性系统在D型学习律下的稳定性和收敛条件作了证明。
2. 初始值问题
运用迭代学习控制技术设计控制器时，只需要通过重复操作获得的受控对象的误差或误差导数信号。在这种控制技术中，迭代学习总要从某初始点开始，初始点指初始状态或初始输出。几乎所有的收敛性证明都要求初始条件是相同的，解决迭代学习控制理论中的初始条件问题一直是人们追求的目标之一。目前已提出的迭代学习控制算法大多数要求被控系统每次运行时的初始状态在期望轨迹对应的初始状态上，即满足初始条件： 

xk(0)=xd(0),k=0,1,2,…（5.24）

当系统的初始状态不在期望轨迹上，而在期望轨迹的某一很小的邻域内时，通常把这类问题归结为学习控制的鲁棒性问题研究。
3. 学习速度问题
在迭代学习算法研究中，其收敛条件基本上都是在学习次数k→∞下给出的。而在实际应用场合，学习次数k→∞显然是没有任何实际意义的。因此，如何使迭代学习过程更快地收敛于期望值是迭代学习控制研究中的另一个重要问题。
ILC本质上是一种前馈控制技术，大部分学习律尽管证明了学习收敛的充分条件，但收敛速度还是很慢。可利用多次学习过程中得到的知识来改进后续学习过程的速度，例如，采用高阶迭代控制算法、带遗忘因子的学习律、利用当前项或反馈配置等方法构造学习律，可使收敛速度大大加快。
4. 鲁棒性问题
迭代学习控制理论的提出有浓厚的工程背景，因此仅仅在无干扰条件下讨论收敛性问题是不够的，还应讨论存在各种干扰的情形下系统的跟踪性能。一个实际运行的迭代学习控制系统除了存在初始偏移外，还或多或少存在状态扰动、测量噪声、输入扰动等各种干扰。鲁棒性问题讨论存在各种干扰时迭代学习控制系统的跟踪性能。具体地说，一个迭代学习控制系统是鲁棒的，指系统在各种有界干扰的影响下，其迭代轨迹能收敛到期望轨迹的邻域内，而当这些干扰消除时，迭代轨迹会收敛到期望轨迹。
5.3机械手轨迹跟踪迭代学习控制仿真实例
5.3.1控制器的设计

考虑一个n关节的机械手，其动态性能可以由以下二阶非线性微分方程描述： 

M(q)q¨+C(q,q·)q·+G(q)=u-τd（5.25）

其中，q∈Rn为关节角位移量，M(q)∈Rn×n为机械手的惯性矩阵，C(q,q·)∈Rn表示离心力和哥氏力，G(q)∈Rn为重力项，τ∈Rn为控制力矩，τd∈Rn为各种误差和扰动。
设系统所要跟踪的期望轨迹为yd(t)，t∈［0,T］。系统第i次输出为yi(t)，令ei(t)=yd(t)-yi(t)。

在学习开始时，系统的初始状态为x0(0)。学习控制的任务为通过学习控制律设计ui+1(t)，使第i+1次运动误差ei+1(t)减少。
采用以下三种基于反馈的迭代学习控制律： 
（1） 闭环D型： 

uk+1(t)=uk(t)+Kd(q·d(t)-q·k+1(t))（5.26）

（2） 闭环PD型： 

uk+1(t)=uk(t)+Kp(qd(t)-qk+1(t))+Kd(q·d(t)-q·k+1(t))（5.27）

（3） 指数变增益D型： 

uk+1(t)=uk(t)+Kp(qd(t)-qk+1(t))+Kd(q·d(t)-q·k+1(t))（5.28）

5.3.2仿真实例
本节针对二关节机械手，介绍一种机械手PD型反馈迭代学习控制的仿真设计方法。针对二关节机械手控制系统式（5.25），各项表示为： 

M= ［mij］2×2

其中，m11=m1l2c1+m2(l21+l2c2+2l1lc2cosq2)+I1+I2，
m12=m21=m2(l2c2+l1lc2cosq2)+l2，
m22=m2l2c2+I2，

C= ［cij］2×2

其中，c11=hq·2，c12=hq·1+hq·2，c21=-hq·1，c22=0，h=-m2l1lc2sinq2

G=［G1G2］T

其中，G1=(m1lc1+m2l1)gcosq1+m2lc2gcos(q1+q2)，G2=m2lc2gcos(q1+q2)，
干扰项为τd=［0.3sint0.1(1-e-t)］T。
机械手系统参数为m1=10，m2=5，l1=1，l2=0.5,lc1=0.5，lc2=0.25,I1=0.83，I2=0.3,g=9.8。
根据式（5.26)~式(5.28），
采用三种闭环迭代学习控制律，其中，M=1为D型迭代学习控制，M=2为PD型迭代学习控制，M=3为指数变增益D型迭代学习控制。
两个关节的角度指令分别为sin(3t)和cos(3t)，为了保证被控对象初始输出与指令初值一致，取被控对象的初始状态为x(0)=［0310］T。取PD型迭代学习控制，即
M=2，仿真结果如图51~图53所示。



图5120次迭代学习的角度跟踪过程





图52第20次迭代学习的角度跟踪




图53第20次迭代学习角度和角速度跟踪误差收敛过程


仿真程序如下： 
(1) 主程序： chap5_1main.m。

clear all;

close all;



t=［0:0.01:3］';

k(1:301)=0;%Total initial points

k=k';

T1(1:301)=0;

T1=T1';

T2=T1;

T=［T1 T2］;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=20;

for i=0:1:M% Start Learning Control

i

pause(0.01);



sim('chap5_1sim',［0,3］);



q1=q(:,1);

dq1=q(:,2);

q2=q(:,3);

dq2=q(:,4);



q1d=qd(:,1);

dq1d=qd(:,2);

q2d=qd(:,3);

dq2d=qd(:,4);



e1=q1d-q1;

e2=q2d-q2;

de1=dq1d-dq1;

de2=dq2d-dq2;



figure(1);

subplot(211);

hold on;

plot(t,q1,'b',t,q1d,'r');

xlabel('time(s)');ylabel('q1d,q1 (rad)');



subplot(212);

hold on;

plot(t,q2,'b',t,q2d,'r');

xlabel('time(s)');ylabel('q2d,q2 (rad)');



j=i+1;

times(j)=i;

e1i(j)=max(abs(e1));

e2i(j)=max(abs(e2));

de1i(j)=max(abs(de1));

de2i(j)=max(abs(de2));

end%End of i

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2);

subplot(211);

plot(t,q1d,'r',t,q1,'b');

xlabel('time(s)');ylabel('angle tracking of link 1');

subplot(212);

plot(t,q2d,'r',t,q2,'b');

xlabel('time(s)');ylabel('angle tracking of link 2');



figure(3);

subplot(211);

plot(t,dq1d,'r',t,dq1,'b');

xlabel('time(s)');ylabel('angle speed tracking of link 1');

subplot(212);

plot(t,dq2d,'r',t,dq2,'b');

xlabel('time(s)');ylabel('angle speed tracking of link 2');



figure(4);

subplot(211);

plot(times,e1i,'*-r',times,e2i,'o-b');

title('change of maximum absolute value of error1 and error2 with times i');

xlabel('time(s)');ylabel('angle tracking error of lmi1 and link 2');

subplot(212);

plot(times,de1i,'*-r',times,de2i,'o-b');

title('change of maximum absolute value of derror1 and derror2 with times i');

xlabel('time(s)');ylabel('angle speed tracking error of lmi1 and link 2');


(2) Simulink子程序： chap5_1sim.mdl。




(3) 被控对象子程序： chap5_1plant.m。

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 1,

sys=mdlDerivatives(t,x,u);

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 4;

sizes.NumDiscStates= 0;

sizes.NumOutputs = 4;

sizes.NumInputs= 2;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0= ［0;3;1;0］;%Must be equal to x(0) of ideal input

str = ［］;

ts= ［0 0］;

function sys=mdlDerivatives(t,x,u)

Tol=［u(1) u(2)］';



g=9.8;

m1=10;m2=5;

l1=1;l2=0.5;

lc1=0.5;lc2=0.25;

I1=0.83;I2=0.3;



d11=m1*lc1^2+m2*(l1^2+lc2^2+2*l1*lc2*cos(x(3)))+I1+I2;

d12=m2*(lc2^2+l1*lc2*cos(x(3)))+I2;

d21=d12;

d22=m2*lc2^2+I2;

D=［d11 d12;d21 d22］;

h=-m2*l1*lc2*sin(x(3));

c11=h*x(4);

c12=h*x(4)+h*x(2);

c21=-h*x(2);

c22=0;

C=［c11 c12;c21 c22］;

g1=(m1*lc1+m2*l1)*g*cos(x(1))+m2*lc2*g*cos(x(1)+x(3));

g2=m2*lc2*g*cos(x(1)+x(3));

G=［g1;g2］;



a=1.0;

d1=a*0.3*sin(t);

d2=a*0.1*(1-exp(-t));

Td=［d1;d2］;



S=-inv(D)*C*［x(2);x(4)］-inv(D)*G+inv(D)*(Tol-Td);



sys(1)=x(2);

sys(2)=S(1);

sys(3)=x(4);

sys(4)=S(2);

function sys=mdlOutputs(t,x,u)

sys(1)=x(1); %Angle1:q1

sys(2)=x(2); %Angle1 speed:dq1

sys(3)=x(3); %Angle2:q2

sys(4)=x(4); %Angle2 speed:dq2
(4) 控制器子程序： chap5_1ctrl.m。

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates= 0;

sizes.NumOutputs = 2;

sizes.NumInputs= 8;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0= ［］;

str = ［］;

ts= ［0 0］;

function sys=mdlOutputs(t,x,u)

q1d=u(1);dq1d=u(2);

q2d=u(3);dq2d=u(4);



q1=u(5);dq1=u(6);

q2=u(7);dq2=u(8);



e1=q1d-q1;

e2=q2d-q2;

e=［e1 e2］';

de1=dq1d-dq1;

de2=dq2d-dq2;

de=［de1 de2］';



Kp=［100 0;0 100］;

Kd=［500 0;0 500］;



M=2;

if M==1

Tol=Kd*de;%D Type

elseif M==2

Tol=Kp*e+Kd*de;%PD Type

elseif M==3

Tol=Kd*exp(0.8*t)*de;%Exponential Gain D Type

end

sys(1)=Tol(1);

sys(2)=Tol(2);
(5) 指令程序： chap5_1input.m。

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates= 0;

sizes.NumOutputs = 4;

sizes.NumInputs= 0;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0= ［］;

str = ［］;

ts= ［0 0］;

function sys=mdlOutputs(t,x,u)

q1d=sin(3*t);

dq1d=3*cos(3*t);

q2d=cos(3*t);

dq2d=-3*sin(3*t);



sys(1)=q1d;

sys(2)=dq1d;

sys(3)=q2d;

sys(4)=dq2d;
5.4线性时变连续系统迭代学习控制
5.4.1系统描述

Arimoto［1］等给出了线性时变连续系统为

x·(t)=A(t)x(t)+B(t)u(t)


y(t)=C(t)x(t)（5.29）

其开环PID型迭代学习控制律为

uk+1(t)=uk(t)+Γddt+L+Ψ∫dtek(t)（5.30）

其中，Γ，L，Ψ为学习增益矩阵。
5.4.2控制器设计及收敛性分析
定理5.1［4］若由式（5.29）和式（5.30）描述的系统满足如下条件： 
（1） ‖I-C(t)B(t)Γ(t)‖≤ρ-<1。
（2） 每次迭代初始条件一致，即xk(0)=x0(k=1,2,3,…),y0(0)=yd(0)，则当k→∞时，有yk(t)→yd(t)，t∈［0,T］。
下面给出该定理的简单分析，可参考文献［4］的证明过程。
由式（5.29）及条件式（5.30）得

yk+1(0)=Cxk+1(0)=Cxk(0)=yk(0)

则ek(0)=0(k=0,1,2,…)，即系统满足初始条件。
非齐次一阶线性微分方程x·(t)=A(t)x(t)+B(t)u(t)的解为

x(t)=Cexp∫t0Adτ+exp∫t0Adτ∫t0B(τ)u(τ)exp∫τ0-Adδdτ

=Cexp(At)+exp(At)∫t0B(τ)u(τ)exp(-Aτ)dτ

=Cexp(At)+∫t0exp(A(t-τ))B(τ)u(τ)dτ

取Φ(t,τ)=exp(A(t-τ))，则
xk(t)-xk+1(t)=∫t0Φ(t,τ)B(τ)(uk(τ)-uk+1(τ))dτ

由于ek(t)=yd(t)-yk(t)，ek+1(t)=yd(t)-yk+1(t)，则

ek+1(t)-ek(t)=yk(t)-yk+1(t)=C(t)(xk(t)-xk+1(t))

=∫t0C(t)Φ(t,τ)B(τ)(uk(τ)-uk+1(τ))dτ

即

ek+1(t)=ek(t)-∫t0C(t)Φ(t,τ)B(τ)(uk+1(τ)-uk(τ))dτ

将PID型控制律式（5.30）代入上式，则第k+1次输出的误差为： 

ek+1(t)=ek(t)-∫t0C(t)Φ(t,τ)B(τ)［Γ(τ)e·k(τ)+L(τ)ek(τ)+Ψ(τ)∫τ0ek(δ)dδ］dτ（5.31）

利用分部积分公式，令G(t,τ)=C(t)B(τ)Γ(τ)，有

∫t0C(t)B(τ)Γ(τ)e·k(τ)dτ=G(t,τ)ek(τ)t0-∫t0τG(t,τ)ek(τ)dτ

=C(t)B(τ)Γ(τ)ek(τ)-∫t0τG(t,τ)ek(τ)dτ（5.32）


将式（5.32）代入式（5.31），得

ek+1(t)=［I-C(t)B(t)Γ(t)］ek(t)+∫t0τG(t,τ)ek(τ)dτ-

∫t0C(t)Φ(t,τ)B(τ)L(τ)ek(τ)dτ-

∫t0∫τ0C(t)Φ(t,τ)B(τ)ψ(τ)ek(σ)dσdτ
（5.33）

将式（5.33）两端取范数，有

‖ek+1(t)‖≤‖I-C(t)B(t)Γ(t)‖‖ek(t)‖+∫t0τG(t,τ)‖ek(τ)‖dτ+

∫t0‖C(t)Φ(t,τ)B(τ)L(τ)‖‖ek(τ)dτ+

∫t0∫τ0C(t)Φ(t,τ)B(τ)ψ(τ)‖‖ek(σ)‖dσdτ

≤‖I-C(t)B(t)Γ(t)‖‖ek(t)‖+∫t0b1‖ek(τ)‖dτ+

∫t0∫τ0b2‖ek(σ)‖dσdτ（5.34）

其中，

b1=maxsupt,τ∈［0,T］τG(t,τ),supt,τ∈［0,T］‖C(t)Φ(t,τ)B(τ)L(τ)‖

b2= supt,τ∈［0,T］‖C(t)Φ(t,τ)B(τ)ψ(τ)‖

将式（5.34）两端同乘以exp(-λt)，λ>0，并考虑∫t0exp(λτ)dτ=exp(λt)-1λ，有

exp(-λt)∫t0b1‖ek(τ)‖dτ=exp(-λt)∫t0b1‖ek(τ)‖exp(-λτ)exp(λτ)dτ

≤b1exp(-λt)‖ek(τ)‖λ∫t0exp(λτ)dτ

=b1exp(-λt)‖ek(τ)‖λexp(λt)-1λ

=b1λ‖ek(τ)‖λexp(-λt)(exp(λt)-1)

=b1(1-exp(-λt))λ‖ek(τ)‖λ

≤b1(1-exp(-λT))λ‖ek(τ)‖λ

根据范数的性质［4］，有

exp(-λt)∫t0∫τ0b2‖ek(σ)‖dσdτ≤b21-exp(-λT)λ2‖ek(τ)‖λ

则

‖ek+1‖λ≤ρ~‖ek‖λ（5.35）

其中，ρ~=ρ-+b11-exp(-λT)λ+b21-exp(-λT)λ2。由于ρ-<1，则当λ取足够大时，可使ρ~<1。因此limk→∞‖ek‖λ=0。
如果将式（5.30）中的e(k)改为e(k+1)，则为闭环PID型迭代学习控制律。同定理5.1的分析过程，可证明闭环PID迭代学习控制律。
5.4.3仿真实例
考虑二输入、二输出线性系统： 

x·1(t)
x·2(t)=-23
11x1(t)
x2(t)+11
01u1(t)
u2(t)

y1(t)
y2(t)=20
01x1(t)
x2(t)

期望跟踪轨迹为
y1d(t)
y2d(t)=sin(3t)
cos(3t)，t∈［0,1］

由于CB=22
01，取Γ=0.950
00.95，可满足定理5.1中的条件（1），在控制律式（5.30）中取L=2.00
02.0，系统的初始状态为x1(0)(0)
x2(0)(0)=0
1。
在chap5_2sim.mdl程序中，选择Simulink的Manual Switch开关，将开关向下，取PD型开环迭代学习控制律，仿真结果如图54~图56所示。将开关向上，采用PD型闭环迭代学习控制律，仿真结果如图57~图59所示。可见，闭环收敛速度好于开环收敛速度。


图5430次迭代学习的跟踪过程（开环PD控制）




图55第30次迭代学习的位置跟踪（开环PD控制）




图5630次迭代学习的跟踪过程（闭环PD控制）




图57第30次迭代学习的位置跟踪（闭环PD控制）




图5830次迭代过程中误差最大绝对值的收敛过程（开环PD控制）





图5930次迭代过程中误差最大绝对值的收敛过程（闭环PD控制）



仿真程序如下： 
（1） 主程序： chap5_2main.m。

%Iterative D-Type Learning Control

clear all;

close all;



t=［0:0.01:1］';

k(1:101)=0;%Total initial points

k=k';

T1(1:101)=0;

T1=T1';

T2=T1;

T=［T1 T2］;



k1(1:101)=0;%Total initial points

k1=k1';

E1(1:101)=0;

E1=E1';

E2=E1;

E3=E1;

E4=E1;

E=［E1 E2 E3 E4］;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=30;

for i=0:1:M% Start Learning Control

i

pause(0.01);



sim('chap5_2sim',［0,1］);



x1=x(:,1);

x2=x(:,2);



x1d=xd(:,1);

x2d=xd(:,2);

dx1d=xd(:,3);

dx2d=xd(:,4);



e1=E(:,1);

e2=E(:,2);

de1=E(:,3);

de2=E(:,4);

e=［e1 e2］';

de=［de1 de2］';



figure(1);

subplot(211);

hold on;

plot(t,x1,'b',t,x1d,'r');

xlabel('time(s)');ylabel('x1d,x1');



subplot(212);

hold on;

plot(t,x2,'b',t,x2d,'r');

xlabel('time(s)');ylabel('x2d,x2');



j=i+1;

times(j)=i;

e1i(j)=max(abs(e1));

e2i(j)=max(abs(e2));

de1i(j)=max(abs(de1));

de2i(j)=max(abs(de2));

end%End of i

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2);

subplot(211);

plot(t,x1d,'r',t,x1,'b');

xlabel('time(s)');ylabel('position tracking of x1');

subplot(212);

plot(t,x2d,'r',t,x2,'b');

xlabel('time(s)');ylabel('position tracking of x2');




figure(3);

plot(times,e1i,'*-r',times,e2i,'o-b');

title('Change of maximum absolute value of error1 and error2 with times');

xlabel('time(s)');ylabel('error1 and error2');

(2) Simulink程序： chap5_2sim.mdl。




(3) 被控对象子程序： chap5_2plant.m。

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 1,

sys=mdlDerivatives(t,x,u);

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 2;

sizes.NumDiscStates= 0;

sizes.NumOutputs = 2;

sizes.NumInputs= 2;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0= ［0;1］;

str = ［］;

ts= ［0 0］;

function sys=mdlDerivatives(t,x,u)

A=［-2 3;1 1］;

C=［1 0;0 1］;

B=［1 1;0 1］;

Gama=0.95;

norm(eye(2)-C*B*Gama);% Must be smaller than 1.0



U=［u(1);u(2)］;

dx=A*x+B*U;

sys(1)=dx(1);

sys(2)=dx(2);

function sys=mdlOutputs(t,x,u)

sys(1)=x(1);

sys(2)=x(2);
(4) 控制器子程序： chap5_2ctrl.m。

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates= 0;

sizes.NumOutputs = 2;

sizes.NumInputs= 4;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0= ［］;

str = ［］;

ts= ［0 0］;

function sys=mdlOutputs(t,x,u)

e1=u(1);e2=u(2);

de1=u(3);de2=u(4);



e=［e1 e2］';

de=［de1 de2］';



Kp=2.0;

Gama=0.95;

Kd=Gama*eye(2); 



Tol=Kp*e+Kd*de;%PD Type



sys(1)=Tol(1);

sys(2)=Tol(2);
(5) 指令程序： chap5_2input.m。

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates= 0;

sizes.NumOutputs = 4;

sizes.NumInputs= 0;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0= ［］;

str = ［］;

ts= ［0 0］;

function sys=mdlOutputs(t,x,u)

x1d=sin(3*t);

dx1d=3*cos(3*t);

x2d=cos(3*t);

dx2d=-3*sin(3*t);



sys(1)=x1d;

sys(2)=x2d;

sys(3)=dx1d;

sys(4)=dx2d;
5.5任意初始状态下的迭代学习控制
下面介绍一种任意初始状态下的学习控制方法［10］及其仿真设计方法。
5.5.1问题的提出
假设一种系统为

x·(t)=Ax(t)+Bu(t)
y(t)=Cx(t)
x(t0)=x(0)(5.36)

其中，x(t)∈Rn，u(t),y(t)∈Rm，A、B、C为相应维数的常阵且满足假设

rank(CB)=m(5.37)

设系统所要跟踪的期望轨迹为yd(t)，t∈［0,T］。系统第i次输出为yi(t)，令ei(t)=yd(t)-yi(t)。
在学习开始时，系统的初始状态为x0(0)，初始控制输入为u0(t)。学习控制的任务为已知第i次运动的ui(t)、xi(0)和ei(t)，通过学习控制律设计ui+1(t)和xi+1(0)，第i+1次运动误差ei+1(t)将减少。
5.5.2控制器的设计
首先介绍范数如下： 

‖e(t)‖∞= max1≤i≤me(i)(t)(5.38)

‖G‖∞= max1≤i≤m∑mj=1g(i,j)(5.39)

‖e(t)‖λ= sup1≤t≤Texp(-λt)‖e(t)‖∞(5.40)

其中，e(i)(t)为e(t)∈Rm中的第i个元素，g(i,j)是G∈Rm×m中的第i,j个元素，λ>0。
学习控制律及初始状态学习律分别为

ui+1(t)=ui(t)+Le·i(t)(5.41)

xi+1(0)=xi(0)+BLei(0)(5.42)

其中，L∈Rm×m为常阵。
定理5.2［10］若学习控制律(5.41)及初始状态学习律(5.42)满足以下条件： 
（1） u0(t)在［0,T］上连续，yd(t)在［0,T］上连续可微。
（2） ‖I∞-CBL‖<1。

则当i→∞时，有
yi(t)→yd(t)，t∈［0,T］
下面给出该定理的详细分析过程，可参考文献［10］的证明过程。
由方程式（5.36）的解为

x(t)=exp(At)x(0)+∫t0exp(A(t-τ))Bu(t)dτ

则

xi+1(t)=exp(At)xi+1(0)+∫t0exp(A(t-τ))Bui+1(τ)dτ

将式（5.41）和（5.42）代入上式，得

xi+1(t)=exp(At)［xi(0)+BLei(0)］+∫t0exp(A(t-τ))［Bui(t)+Le·i(t)］dτ

则

ei+1(t)=yd(t)-yi+1(t)=yd(t)-Cxi+1(t)

=yd(t)-Cexp(At)［xi(0)+BLei(0)］+∫t0exp(A(t-τ))［Bui(t)+Le·i(t)］dτ

=yd(t)-Cexp(At)xi(0)+Cexp(At)BLei(0)+

∫t0Cexp(A(t-τ))Bui(t)dτ+∫t0Cexp(A(t-τ))BLe·i(t)dτ

采用分部积分方法： 

∫t0xy·dt=xy|t0-∫t0x·ydt

则

∫t0Cexp(A(t-τ))BLe·i(t)dτ

=［Cexp(A(t-τ))BLei(t)］t0-∫t0(-1)CAexp(A(t-τ))BLei(t)dτ

=CBLei(t)-Cexp(At)BLei(0)+∫t0CAexp(A(t-τ))BLei(τ)dτ

ei+1(t)=yd(t)-Cexp(At)xi(0)-Cexp(At)BLei(0)-∫t0Cexp(A(t-τ))Bui(τ)dτ-

CBLei(t)+Cexp(At)BLei(0)-∫t0CAexp(A(t-τ))BLei(τ)dτ

=yd(t)-Cexp(At)xi(0)+∫t0exp(A(t-τ))Bui(τ)dτ-

CBLei(t)-∫t0CAexp(A(t-τ))BLei(τ)dτ

=yd(t)-yi(t)-CBLei(t)-∫t0CAexp(A(t-τ))BLei(τ)dτ

即

ei+1(t)=ei(t)-CBLei(t)-∫t0CAexp(A(t-τ))BLei(τ)dτ

=(Im-CBL)ei(t)-∫t0CAexp(A(t-τ))BLei(τ)dτ

上式两边同乘以e-λt，取范数，并考虑
‖XY‖≤‖X‖‖Y‖，‖X+Y‖≤‖X‖+‖Y‖

则

 exp(-λt)‖ei+1(t)‖∞

≤‖(Im-CBL)ei(t)‖∞exp(-λt)+∫t0CAexp(A(t-τ))BLei(τ)dτ∞exp(-λt)

≤‖(Im-CBL)‖∞‖ei(t)‖∞exp(-λt)+∫t0CAexp(A(t-τ))BLei(τ)dτ∞exp(-λt)

=‖(Im-CBL)‖∞‖ei(t)‖λ+‖CABL‖∞∫t0exp(A(t-τ))ei(τ)dτ∞exp(-λt)

=‖(Im-CBL)‖∞‖ei(t)‖λ+‖CABL‖∞‖h(t)‖λ

其中，h(t)=∫t0exp(A(t-τ))ei(τ)dτ。
根据λ范数的性质［4］，当λ>A时，有

‖h(t)‖λ≤1-exp((A-λ)T)λ-A‖ei(t)‖λ

则

‖ei+1(t)‖λ≤‖Im-CBL‖∞+‖CABL‖∞1-exp((A-λ)T)λ-A‖ei(t)‖λ=ρ‖ei(t)‖λ

定义

ρ=‖Im-CBL‖∞+‖CABL‖∞1-exp((A-λ)T)λ-A(5.43)

当λ足够大时，1-exp((A-λ)T)λ-A→0，考虑条件（2），则
ρ<1
当i→∞时，有
‖ei+1(t)‖λ→0，t∈［0,T］
由ρ的定义可知，当取L=(CB)-1时，‖Im-CBL‖∞=0，ρ的值最小，收敛速度最快。
5.5.3仿真实例
考虑多输入多输出非线性系统： 

x·1(t)
x·2(t)=-23
11x1(t)
x2(t)+11
01u1(t)
u2(t)

y1(t)
y2(t)=20
01x1(t)
x2(t)

期望跟踪轨迹为
y1d(t)
y2d(t)=1.5t
1.5t，t∈［0,1］

由于CB=22
01，故取L=(CB)-1=0.5-1
01，可满足定理5.2中的条件（2）。
根据式（5.41）和式（5.42），
学习控制律及初始状态学习律分别为

u1(i+1)(t)
u2(i+1)(t)=u1(i)(t)
u2(i)(t)+0.4-0.5
00.5e·1(i)(t)
e·2(i)(t)

x1(i+1)(0)
x2(i+1)(0)=x1(i)(0)
x2(i)(0)+0.40
00.5e1(i)(0)
e2(i)(0)

系统的初始控制输入为u1(0)(t)
u2(0)(t)=0
0，系统的初始状态为x1(0)(0)
x2(0)(0)=2
1。仿真结果如图510~图515所示。


图5105次迭代对象输出的跟踪过程




图5115次迭代后正弦位置跟踪




图5125次迭代过程中误差范数的收敛过程




图51310次迭代对象输出的跟踪过程





图51410次迭代后正弦位置跟踪




图51510次迭代过程中误差范数的收敛过程


仿真程序如下： 
（1） 主程序： chap5_3.m。

%Learning Control with an arbitrary initial state

clear all;

close all;

global A B



A=［-2 3;1 1］;

B=［1 1;0 1］;

C=［2 0;0 1］;

L=inv(C*B);



ts=0.01;

for k=1:1:101

u1(k)=0;u2(k)=0;



end

xk0=［2;1］;

%xk0=［-2;-1］;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=10;

for i=0:1:M% Start Learning Control

i

pause(0.005);

if i==0

xki=xk0;

else

yd0=0;

yi0=［2*xk0(1);xk0(2)］;

e0=yd0-yi0;

xki=xk0+B*L*e0;

end

xk0=xki;%用xk0存储上次运行的初始状态



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:1:101

time(k)=(k-1)*ts;



S=2;

if S==1

y1d_1(k)=1.5*(k-1)*ts;

y2d_1(k)=1.5*(k-1)*ts;

y1d(k)=1.5*k*ts;

y2d(k)=1.5*k*ts;

elseif S==2

y1d_1(k)=sin(4*pi*(k-1)*ts);

y2d_1(k)=sin(4*pi*(k-1)*ts);

y1d(k)=sin(4*pi*k*ts);

y2d(k)=sin(4*pi*k*ts);

end



TimeSet=［(k-1)*ts k*ts］;

para=［u1(k);u2(k)］;

if k==1 % Initial state at times M

xk=xk0;

end



y1_1(k)=2*xk(1);

y2_1(k)=xk(2);

%xk

［tt,xx］=ode45('chap5_3plant',TimeSet,xk,［］,para);

%xx

xk=xx(length(xx),:);

y1(k)=2*xk(1);

y2(k)=xk(2);



e1_1(k)=y1d_1(k)-y1_1(k);

e2_1(k)=y2d_1(k)-y2_1(k);



e1(k)=y1d(k)-y1(k);

e2(k)=y2d(k)-y2(k);



de1(k)=(e1(k)-e1_1(k))/ts;

de2(k)=(e2(k)-e2_1(k))/ts;

dek=［de1(k);de2(k)］; 

Uk=［u1(k);u2(k)］;

U=Uk+L*dek;% Control law: Uk is U(i-1), dek is near to de(i-1)



u1(k)=U(1);

u2(k)=U(2);

end% End of k



figure(1);

subplot(211);

hold on;

plot(time,y1d_1,'r',time,y1_1,'b');

xlabel('time(s)');ylabel('y1d,y1');



subplot(212);

hold on;

plot(time,y2d_1,'r',time,y2_1,'b');

xlabel('time(s)');ylabel('y2d,y2');



i=i+1;

times(i)=i-1;

e1i(i)=max(abs(e1_1));

e2i(i)=max(abs(e2_1));

end%End of i

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2);

subplot(211);

plot(time,y1d_1,'r',time,y1_1,'b');

xlabel('time(s)');ylabel('y1d,y1');

subplot(212);

plot(time,y2d_1,'r',time,y2_1,'b');

xlabel('time(s)');ylabel('y2d,y2');



figure(3);

plot(times,e1i,'*-r',times,e2i,'o-b');

title('change of maximum absolute value of error1 and error2 with times i');

xlabel('times');ylabel('error1 and error2');
(2) 子程序： chap5_3plant.m。

function dx=PlantModel(t,x,flag,para)

global A B

dx=zeros(2,1);



u=para;

dx=A*x+B*u;
实际工程中，机械手常在有限时间内执行重复性的控制任务。针对这种有限区间执行重复性的控制任务，迭代学习控制是一种有效的控制方法。因此，机械臂轨迹跟踪的迭代学习控制得到广泛研究［11~15］。
常规的机械手迭代学习控制方法需要假设机械臂初始状态与期望轨迹初始状态相等，而这在实际过程中常常无法满足，因而针对带有初始角度偏移的机械手迭代学习控制问题的研究具有一定工程意义［14,15］。
参考文献


［1］Arimoto S,Kawamura S,Miyazaki F.Bettering operation of robotics by leaning［J］.Journal of Robotic System,1984,1(2)： 123140.

［2］林辉，王林.迭代学习控制理论［M］.西安： 西北工业大学出版社,1998.

［3］孙明轩，黄宝健.迭代学习控制［M］.北京： 国防工业出版社,1999.

［4］谢胜利.迭代学习控制的理论与应用［M］.北京： 科学出版社，2005.

［5］于少娟，齐向东，吴聚华.迭代学习控制理论及应用［M］.北京： 机械工业出版社，2005.

［6］方忠，韩正之，陈彭年.迭代学习控制新进展［J］.控制理论与应用，2002，19(2)： 161165.

［7］石成英，林辉.迭代学习控制技术的原理、算法及应用［J］.机床与液压，2004，19： 8083.

［8］许建新，侯忠生.学习控制的现状与展望［J］.自动化学报，2005，31(6)： 943955.

［9］李仁俊，韩正之.迭代学习控制综述［J］.控制与决策，2005，20（9）： 961966.

［10］任雪梅，高为炳.任意初始状态下的迭代学习控制［J］.自动化学报，1994，20(1)： 7479.

［11］Ouyang P R,Zhang W J,Gupta M M.An adaptive switching learning control method for trajectory tracking of robot manipulators［J］.Mechatronics,2006,16： 5161.

［12］Tayebi A.Adaptive iterative learning control for robot manipulators［J］.Automatica,2004,40： 11951203.

［13］Kang M K,Lee J S,Han K L.Kinematic pathtracking of mobile robot using iterative learning control［J］.Journal of Robotic Systems,2005,22(2)： 111121.


［14］Chen Y Q,Wen C,Gong Z,et al.An iterative learning controller with initial state learning［J］.IEEE Transaction on Automatic Control,1999,44(2)： 371376.


［15］Park K H,Bien Z.A generalized iterative learning controller against initial state error［J］.International Journal of Control，2000,73(10)： 871881.