随机控制系统仿真

随机控制系统仿真按照原则，控制系统可分为确定性系统和随机系统两大类。确定性系统：系统参数、指令、干扰均为已知量随机系统：系统参数、指令、干扰都是随机的随机系统在一定条件下可以简化为确定性系统。如果系统的一次行为不能代表多次运行性能，则该系统不能简化为确定性系统，此时，必须按随机系统进行仿真。多数武器系统为随机系统。引言随机系统仿真的目的是获取状态变量、输出等的统计特性。第2页,共65页，2024年2月25日，星期天7.1概率和随机过程一、随机控制系统系统模型：若a、b为零阵，x0为已知确定值，则该系统为确定性系统。若a、b、u（t）、x0均为不确定的，则该系统为随机系统。第3页,共65页，2024年2月25日，星期天二、基本概念1、随机事件：可能发生，可能不发生的事件。2、随机过程：随机系统状态变量和输出的变化过程。这是一组曲线。既含有一次实验中x随t变化的过程，也包含多次实验在相同时刻x的值。3、样本函数：随机过程任一条状态变量或输出变量曲线。

{xi(t)}Li=1是随机过程。

xk(t)、x1(t)是样本函数4、随机变量：某时刻的状态变量或输出变量称为随机变量。第4页,共65页，2024年2月25日，星期天随机过程样本函数随机变量第5页,共65页，2024年2月25日，星期天三、连续随机变量的统计特性1、概率密度p(y)：每个y值发生的可能性。2、概率分布函数：变量值小于或等于y的随机变量的概率。3、均值（一阶矩/数学期望）：E（y）E(y1+y2+…+yn)=E(y1)+E(y2)+…+E(yn)第6页,共65页，2024年2月25日，星期天4、均方值E（y2）与均方根sqrt（E(y2)）5、方差σ2y与方差根σy：表示随机变量偏离均值的程度。方差又称二阶中心矩。

σ2y=E(y-E(y))2=E(y2)-E2(y)σ2(y1+y2+…+yn)=σ2(y1)+σ2(y2)+…+σ2(yn)随机系统仿真的目的是获取状态变量、输出等的统计特性。第7页,共65页，2024年2月25日，星期天

用仿真方法获取随机变量的统计特性时，随机变量只有有限个数值，此时，各统计特性的计算公式为：均值：均方值：方差：第8页,共65页，2024年2月25日，星期天四、典型概率分布1、均匀分布概率密度：

p(y)=1/(b-a)，a<=y<=b=0，其他概率分布函数：

P(y)=0，y<=a=(y-a)/(b-a)，a<y<=b=1，y>b

均值：E（y）=(a+b)/2，方差：σ2y=(b-a)2/12。第9页,共65页，2024年2月25日，星期天2、正态分布：自然界中最常见的一种分布。又称高斯分布，钟形分布。σ2y为方差，m为均值。概率密度：第10页,共65页，2024年2月25日，星期天中心极限定理：正态概率分布随机变量可以用无数个任意分布的相互独立的随机变量合成。在工程上，可用10~12个互相独立的均匀分布合成一个正态分布。概率分布函数：第11页,共65页，2024年2月25日，星期天五、随机过程的统计特性平稳随机过程：统计特性不随时间变化。各态历经平稳随机过程：一个样本函数的时间统计特性等于随机过程统计特性。常见的随机过程都是各态历经的平稳随机过程，故可用样本函数的统计特性来表征随机过程的统计特性。第12页,共65页，2024年2月25日，星期天1、均值：yk(t)为样本，Y(t)为随机过程。2、方差：第13页,共65页，2024年2月25日，星期天3、自相关函数：衡量随机过程功率强弱的尺度。反映Y（t1）与Y（t2）的相关性。RY(t1，t2)大，则Y(t)变化平缓，可预见性大，RY(t1，t2)小，则Y(t)变化剧烈，可预见性小。第14页,共65页，2024年2月25日，星期天4、功率谱密度（能量谱密度）是RY(t1，t2)的傅氏变换，反映Y（t）中不同频率的能量。5、白色随机过程：

最简单的随机过程

最理想的随机过程均值：E（yk）=0，方差：σ2，自相关函数：Ry(τ)=σ2δ(τ)

功率谱密度：Sy(ω)=σ2，各频率分量能量相等，带宽无限。白色随机过程的自相关函数及功率谱密度第15页,共65页，2024年2月25日，星期天六、随机过程作用下线性系统的响应1、输入输出的时域关系系统模型为脉冲响应函数h（t）y(t)与x(t)为卷积关系均方值：表明输出的均方值与输入的自相关函数有关对白色随机过程：第16页,共65页，2024年2月25日，星期天2、输入输出的频域关系系统模型为传递函数Sy(ω)=|H(jω)|2

Sx(ω)若Sx(ω)=σ2

，|H(jω)|2不为常值，则y(t)为有色噪声。例：H(jω)=1/（jωT+1)，Sx(ω)=σ2

则，y(t)为有色噪声。第17页,共65页，2024年2月25日，星期天白色随机过程作用下线性系统的频域响应第18页,共65页，2024年2月25日，星期天7.2随机控制系统仿真的专门问题随机数产生随机数概率分布的转变随机数概率功率谱密度的转变系统随机参数误差的产生系统随机初值的产生系统随机干扰的产生第19页,共65页，2024年2月25日，星期天一、随机数产生产生伪随机数序列=>改变概率分布、功率谱密度=>加入随机系统=>进行仿真1、公式：计算机产生伪随机序列，乘同余法公式：x(i+1)=(λ*x(i))(modμ)其中：λ=2k/2，且λ=8N+/-3，K为字长，N为整数。

μ=2n，且0<μ<2k，n为整数2、特点（1）伪随机序列，周期为μ/4（2）均匀分布，值域：0—2n-1第20页,共65页，2024年2月25日，星期天3、统计特性（1）概率密度：p（x）=1/μ，0<=x<=μ=0，x〈0或x〉μ（2）均值：E（x）=μ/2（3）方差：σ2x=μ2/12（4）自相关函数：Rx(τ)=0，|τ|>T=μ2（1-τ/T）/12，|τ|<=T第21页,共65页，2024年2月25日，星期天当ω较小时，Sx（ω

）近似为常值（5）功率谱密度：在ω<π/T频带内，可将随机数序列看作白色随机过程。第22页,共65页，2024年2月25日，星期天例：设计算机的字长为32位，用乘同余法产生随机数序列。解：取μ=2^30，N=8192，λ=8*N+3=65539，x（1）=11111。

编程序计算随机数序列。运行该程序，得到随机数序列。随机数序列随机数序列（局部放大）第23页,共65页，2024年2月25日，星期天令随机数产生周期为T=0.01s，画出其时间曲线，如下图。随机数的时间曲线随机数的时间曲线（局部）第24页,共65页，2024年2月25日，星期天乘同余法产生的伪随机序列的自相关函数由图可知，

τ=0时，Rx(τ)最大，为9.6077e+016，即：μ2T/12τ>T时，Rx(τ)=0其自相关函数的最大值为：Rx（0）=μ2T/12=9.6077e+016τ>0.01时，Rx(τ)=0第25页,共65页，2024年2月25日，星期天乘同余法产生的伪随机序列的功率谱密度功率谱密度：Sx（0）=μ

^2*T/12=9.607679205057059e+014低频时，Sx（ω

）近似为常值。T越小，为常值的频带越大。第26页,共65页，2024年2月25日，星期天二、随机数概率分布的改变多数情况下需要将均匀分布改为正态分布。中心极限定理：正态概率分布随机变量可以用无数个任意分布的相互独立的随机变量合成。在工程上，可用10~12个互相独立的均匀分布合成一个正态分布。第27页,共65页，2024年2月25日，星期天1、产生值域为（0，2k-1）的均匀分布随机数序列{x1(i)}。2、将该序列标称化，值域变为（0，1）方法：x1(i)*1/(2k-1)3、均值移位，值域变为（-0.5，0.5）方法：x1(i)-0.5此时，E（x）=0，σ2x=1/12

重复1-3步，产生n个不相关的均匀分布随机数序列，记为{x2(i)}、{x3(i)}……{xn(i)}

。4、用上述n个不相关的均匀分布随机数序列构成一个正态分布随机数序列。第28页,共65页，2024年2月25日，星期天序列1：x11,x21,……xm1序列2：x12,x22,……xm2...序列n：x1n,x2n,……xmn一般，取n=12，则合成的随机数序列E（x）=0，σ2x=1……合成：x1,x2,……xn第29页,共65页，2024年2月25日，星期天三、随机数功率谱密度的改变利用线性系统做成形滤波器Sy(ω)=|H(jω)|2*μ2T/12，其中，Sy(ω)为需要的功率谱密度，

μ2T/12为伪随机序列的功率谱密度。

则根据Sy(ω)可求出线性系统的频率特性H(jω)。第30页,共65页，2024年2月25日，星期天四、仿真中一些随机参数的产生1、参数误差及随机初值的产生

二者均为正态分布，所给参数为最大允许误差ΔA。任务：根据最大允许误差产生随机参数误差或随机初值序列。第31页,共65页，2024年2月25日，星期天分析：由正态分布的性质，一般认为最大允许误差

ΔA>=4σ，

σ为随机序列的方差根。即：随机序列的方差为（取等号情况）：

σ2=（ΔA/4）2（1）若产生12个值域为（-b，b）的独立的均匀分布随机序列，其方差均为

σ2=b2/3则合成后所得正态分布序列的方差为

σ2=4b2（2)考虑（1）式，则有

b=ΔA/8即：所产生的均匀分布序列的值域为（-ΔA/8，ΔA/8）第32页,共65页，2024年2月25日，星期天2、随机干扰的产生随机干扰既要求概率统计特性，又要求频谱特性。一般认为随机干扰的频谱可通过实测或分析得到，即是已知的，故可通过下述方法产生随机干扰。第一步，产生白色随机过程（在一定频带内）。第二步，根据随机干扰的频谱设计成形滤波器。第三步，将白色随机过程作用于成形滤波器，得到有色随机过程。第33页,共65页，2024年2月25日，星期天7.3MonteCarlo法

适用于线性或非线性随机系统仿真，使用时限制条件少，但仿真工作量大。通过对仿真所得数据进行处理得到系统状态变量和输出的统计特性。第34页,共65页，2024年2月25日，星期天蒙特卡罗方法的基本思想

二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，人们在生产实践和科学试验中就已发现，并加以利用。两个例子

例1.蒲丰氏问题

例2.射击问题（打靶游戏）基本思想第35页,共65页，2024年2月25日，星期天例1.蒲丰氏问题

为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（

l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：

求出π值其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。第36页,共65页，2024年2月25日，星期天

一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929第37页,共65页，2024年2月25日，星期天例2.射击问题（打靶游戏）

设r表示射击运动员的弹着点到靶心的距离，ｇ(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为

用概率语言来说，<g>是随机变量ｇ(r)的数学期望，即第38页,共65页，2024年2月25日，星期天

现假设该运动员进行了Ｎ次射击，每次射击的弹着点依次为r1，r2，…，rN，则Ｎ次得分g(r1)，g(r2)，…，g(rN)的算术平均值代表了该运动员的成绩。换言之，为积分<g>的估计值，或近似值。在该例中，用Ｎ次试验所得成绩的算术平均值作为数学期望<g>的估计值（积分近似值）。

第39页,共65页，2024年2月25日，星期天基本思想

由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。第40页,共65页，2024年2月25日，星期天

因此，可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量ｇ(r)的数学期望

通过某种试验，得到Ｎ个观察值r1，r2，…，rN（用概率语言来说，从分布密度函数f(r)中抽取Ｎ个子样r1，r2，…，rN，），将相应的Ｎ个随机变量的值g(r1)，g(r2)，…，g(rN)的算术平均值

作为积分的估计值（近似值）。第41页,共65页，2024年2月25日，星期天

为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。

第42页,共65页，2024年2月25日，星期天一、蒙特卡罗法仿真步骤1、建立随机系统模型，产生随机参数序列。2、改变随机参数的值进行仿真。多采用数值积分法。一般，每个随机参数都要取不同的值仿真50次以上，以使结果可信。3、处理仿真结果，求状态变量、输出的统计特性。第43页,共65页，2024年2月25日，星期天二、仿真结果处理设对某参数共仿真N次，每次计算M步。{yi1}：y11,y21,……yM1{yi2}：y12,y22,……yM2

...{yiN}：y1n,y2n,……yMN第44页,共65页，2024年2月25日，星期天1、均值：2、方差：根据仿真结果求统计特性：均值、方差、概率密度及频谱。第45页,共65页，2024年2月25日，星期天3、概率密度：设Δyi为{yi}值域的分组区间，Δyi=(b–a)/LNyij为介于区间(yij-0.5Δyi，yij+0.5Δyi)中yij的个数，N为{yi}中yij的总数，即仿真次数。则，概率密度为：pi(yij)=Nyij/NΔyi4、频谱：其中：j为虚单位，r为谐波次数，Ar为幅值。第46页,共65页，2024年2月25日，星期天例

dotX(t)=(A+a)X(t)+(B+b)u(t)，

y(t)=CX(t)初值：X0+ΔX0ΔX0为初值误差，正态分布，最大误差为ΔA用MonteCarlo法仿真，求初值误差ΔX0对输出y(t)统计特性的影响。步骤1：产生正态分布的随机数序列ΔX0；用乘同余法产生12个值域为（-ΔA/8，ΔA/8）的均匀分布序列，然后合成即为正态分布的随机数序列ΔX0。步骤2：取不同的ΔX0进行仿真。采用数值积分法，仿真次数在50次以上。步骤3：数据处理，得出初值误差对输出的影响。第47页,共65页，2024年2月25日，星期天仿真流程图第48页,共65页，2024年2月25日，星期天蒙特卡罗法的特点优点(1)能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。(2)受几何条件限制小。(3)收敛速度与问题的维数无关。(4)具有同时计算多个方案与多个未知量的能力。(5)误差容易确定。(6)程序结构简单易于实现。缺点(1)收敛速度慢。(2)误差具有概率性。第49页,共65页，2024年2月25日，星期天能够比较逼真地描述具有随机性质的事物的特点及物理实验过程

从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。第50页,共65页，2024年2月25日，星期天受几何条件限制小

在计算s维空间中的任一区域Ds上的积分无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个，得到积分的近似值。其中Ds为区域Ds的体积。这是数值方法难以实现的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。第51页,共65页，2024年2月25日，星期天收敛速度与问题的维数无关

由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。第52页,共65页，2024年2月25日，星期天具有同时计算多个方案与多个未知量的能力（1）对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。（2）另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。

第53页,共65页，2024年2月25日，星期天误差容易确定

（1）对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。（2）一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。第54页,共65页，2024年2月25日，星期天程序结构简单，易于实现

在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。第55页,共65页，2024年2月25日，星期天7.4伴随系统仿真法一、伴随系统设系统的状态方程为dotX(t)=A(t)X(t)+B(t)u(t)y(t)=C(t)X(t)X(0)=X0，0<=t<=tf，

A：nxn，B：nxm，C：rxn则其伴随系统为：第56页,共65页，2024年2月25日，星期天：伴随状态，nx1：伴随输入，rx1：伴随输出，mx1：伴随系统时间自变量第57页,共65页，2024年2月25日，星期天二、伴随系统结构图