第3 章 科学运算问题的MATLAB求解 　　控制系统的研究涉及各种各样的数学问题求解。例如，系统的稳 定性分析需要求取矩阵的特征根，可控性与可观测性的判定需要求出 矩阵的秩，而状态转移矩阵的求解需要矩阵指数的求解，这需要线性 代数问题的求解。控制系统的仿真需要微分方程的求解，最优控制器 设计涉及最优化问题的求解。如果能掌握利用MATLAB 求解数学问 题的方法和技术无疑会提高控制系统分析与设计的能力。 　　本章将介绍MATLAB 语言在现代科学运算领域中的应用。这里 所谓“运算”是指利用MATLAB 不但能进行数值计算，还可以进行解 析运算，所以从涵盖范围上比一般数值“计算”更广泛。MATLAB 起源 于线性代数的数值运算，在其长期的发展过程中，形成了几乎所有应 用数学分支的求解函数与专门工具箱，并成功地引入了符号运算的功 能，使得公式推导成为可能。所以一般情况下用一个语句就能直接获 得数学问题的解。 　　MATLAB 语言求解科学运算的功能是其广受科学工作者喜爱的 重要原因，也是MATLAB 语言的一大重要的特色。本章侧重于介绍用 MATLAB 为主要工具，直接求解和控制学科密切相关的数学问题的 方法，为更好地探讨控制问题打下良好的基础。 　　3.1 节介绍线性代数问题的解析与数值求解方法，介绍包括矩阵 基本分析、矩阵变换的解析与数值解方法，并介绍矩阵函数的计算。 3.2 节介绍各种方程求解方法，包括线性代数方程、非线性方程和矩 阵方程的求解方法，并试图求出多解方程全部的根。3.3 节介绍一阶 显式微分方程组的数值解方法，并介绍将一般常微分方程变换成可求 解标准型的方法，还将介绍一般线性常系数微分方程的解析解方法。 3.4 节介绍最优化问题的数值求解方法，包括无约束最优化问题、有约 束最优化问题的求解方法和曲线的最小二乘拟合方法。3.5 节将介绍 Laplace 变换与z 变换问题MATLAB 求解方法。还将通过Laplace 变 换数值求解方法得出无理闭环系统的时域响应数值解。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 本章涉及了大量数学公式，但核心问题是引导读者如何避开数学问题本身及 烦琐的底层解法，在MATLAB 框架下直接获得可靠的解。关于应用MATLAB 求 解各种各样数学问题的详细内容可以参阅文献[1～5]。 3.1 线性代数问题的MATLAB求解 线性代数与矩阵运算是控制领域最重要的数学基础之一。很多线性代数问题 是可以求取解析解的，不能求取解析解的问题往往也能得出数值解。本节将以数值 解的介绍为主，其中很多函数同样可以利用MATLAB 的符号运算工具箱中提供的 相同函数名求出解析解。 3.1.1 矩阵的基本分析 矩阵的基本分析往往可以反映出矩阵的某些性质，比如在控制系统分析中，矩 阵的特征值可以用来分析系统的稳定性，矩阵的秩可以用来分析系统的可控性和 可观测性等，这里将系统地介绍矩阵基本分析的概念及其MATLAB 实现。 1. 矩阵的行列式（determinant） MATLAB 提供了内核函数det(A)，利用它可以直接求取矩阵A的行列式。 如果矩阵A为数值矩阵，则得出的行列式为数值计算结果；如果A定义为符号矩 阵，则det() 函数将得出解析解。二者的区别是，对接近奇异的系统来说，解析解方 法得出的结果更精确。不过，解析解算法也有一定的局限性，不适合大规模矩阵的 求解。文献[3] 通过实验探讨了随机符号矩阵的行列式求解，指出如果矩阵不含有 变量，则可以在1 s 内获得30 × 30 符号矩阵行列式的解析解，而90 × 90 矩阵耗时接 近一分钟。 例3-1 试求Hilbert 矩阵的行列式。 解Hilbert 矩阵的通项为hi,j = 1/(i + j ? 1)，用MATLAB 的命令hilb(n) 函数就可 以在MATLAB 工作空间中定义出来，而用sym() 函数即可得出其符号型表示。下面的 语句即可生成并计算出10 阶Hilbert 矩阵的行列式为： >> H=hilb(10); d1=det(H) % 求解数值解 H=sym(H); d2=det(H) % 先将矩阵变换成符号矩阵，再求解析解 用第1 行语句可以求出数值解为d1 = 2.1644×10?53，该结果不精确，故需要用解析解 方法求解。由第2 行命令可以得出解析解为 d2 = 1 46206893947914691316295628839036278726983680000000000 例3-2 试求下面带有变量的Vandermonde 矩阵的行列式。 A = 24 1 1 1 a b c a2 b2 c2 35 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解75 解先定义符号变量a、b、c，然后用下面的语句输入矩阵，并得出矩阵的特征多项式。 >> syms a b c % syms 命令可以声明符号变量，用空格分隔 A=[1,1,1; a,b,c; a^2,b^2,c^2]; % 建立Vandermonde 矩阵 det(A); simplify(factor(ans)) % 求行列式并进行因式分解 可以得出行列式的因式分解形式为?(a ? b)(a ? c)(b ? c)。 2. 矩阵的迹（trace） 假设一个方阵为A = {aij}，则矩阵A的迹定义为该矩阵对角线上各个元素之 和。由代数理论可知，矩阵的迹和该矩阵的特征值之和是相同的，矩阵A的迹可以 由MATLAB 函数trace(A) 求出。 3. 矩阵的秩（rank） 若矩阵所有的列向量中最多有rc 列线性无关，则称矩阵的列秩为rc，如果rc = m，则称A为列满秩矩阵。相应地，若矩阵A的行向量中有rr 个是线性无关的，则 称矩阵A的行秩为rr。如果rr = n，则称A为行满秩矩阵。可以证明，矩阵的行 秩和列秩是相等的，故称为矩阵的秩，记作rank(A) = rc = rr，这时矩阵的秩 为rank(A)。矩阵的秩也表示该矩阵中行列式不等于0 的子式的最大阶次，所谓子 式，即为从原矩阵中任取k 行及k 列所构成的子矩阵。MATLAB 提供了一个内核函 数rank(A,ε) 来用数值方法求取一个已知矩阵A的数值秩，其中ε 为机器精度。如 果没有特殊说明，可以由rank(A) 函数求出A矩阵的秩。 4. 矩阵的范数（norm） 矩阵的常用范数定义为 ||A||1= max 1?j?n n Xi=1 |aij |, ||A||2=psmax(ATA), ||A||∞= max 1?i?n n Xj=1 |aij | （3-1-1） 其中，s(X) 为X矩阵的特征值，而smax(ATA) 即为ATA矩阵的最大特征值。事实 上，||A||2 为A矩阵的最大奇异值。MATLAB 提供了求取矩阵范数的函数norm(A) 可以求出||A||2，矩阵的1-范数||A||1 可以由norm(A,1) 求解，矩阵无穷范数||A||∞ 可以由norm(A,inf) 求出。注意，该函数只能求数值解。 5. 矩阵的特征多项式、特征方程与特征根（eigenvalues） 矩阵sI ? A的行列式可以写成一个关于s 的多项式C(s) C(s) = det(sI ? A) = sn + c1sn?1 + · · · + cn?1s + cn （3-1-2） 这样的多项式C(s) 称为矩阵A的特征多项式，其中系数ci, i = 1, 2, · · · , n称为矩 阵的特征多项式系数。 MATLAB 提供了求取矩阵特征多项式系数的函数p=poly(A)，而返回的p 为一个行向量，其各个分量为矩阵A的降幂排列的特征多项式系数。该函数的另外 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 一种调用格式是：如果给定的A为向量，则假定该向量是一个矩阵的特征根，由此 求出该矩阵的特征多项式系数，如果向量A中有无穷大或NaN 值，则首先剔除。 对方阵A，如果存在一个非零的向量x，且有一个标量λ满足Ax=λx，则称λ 为A矩阵的一个特征值，而x为对应于特征值λ的特征向量，严格说来，x应该称为 A的右特征向量。如果矩阵A的特征值不包含重复的值，则对应的各个特征向量为 线性无关的，这样由各个特征向量可以构成一个非奇异的矩阵，如果用它对原始矩 阵作相似变换，则可以得出一个对角矩阵。矩阵的特征值与特征向量由MATLAB 提供的函数eig() 可以容易地求出，该函数的调用格式为[V ,D]=eig(A)，其中 A为给定的矩阵，解出的D为一个对角矩阵，其对角线上的元素为矩阵A的特征 值，而每个特征值对应于V 矩阵中的一列，称为该特征值的特征向量。MATLAB 的 矩阵特征值的结果满足AV = V D，且V 矩阵每个特征向量各元素的平方和（即 列向量的2 范数）均为1。如果调用该函数时至多只给出一个返回变元，则将返回矩 阵A的特征值。即使A为复数矩阵，也照样可以由eig() 函数得出其特征值与特征 向量矩阵。 6. 多项式及多项式矩阵的求值 可以由C=polyval(a,x) 命令求取基于点运算的多项式值，求出C=a1x.^n +a2x.^(n ? 1) + · · · + an+1，其中a=[a1,a2,· · · ,an,an+1] 为多项式系数降幂排 列构成的向量，x为一个任意的向量或矩阵。 如果想求取真正的矩阵多项式的值，亦即 B = a1An + a2An?1 + · · · + anA + an+1I （3-1-3） 其中，I 为和A同阶次的单位矩阵，则可以用B=polyvalm(a,A)。 7. 矩阵的逆与广义逆 对一个已知的n × n非奇异方阵A，如果有一个C 矩阵满足 AC = CA = I （3-1-4） 其中I 为单位矩阵，则称C 矩阵为A矩阵的逆矩阵，并记作C = A?1。MATLAB 下提供的C=inv(A) 函数即可求出矩阵A的逆矩阵C。 如果用户想得出奇异矩阵或长方形矩阵的一种“逆”阵，则需要使用广义逆的 概念。对一个给定的矩阵A，存在一个唯一的矩阵M同时满足3 个条件： （1）AMA = A。 （2）MAM =M。 （3）AM 与MA均为对称矩阵。 这样的矩阵M称为矩阵A的Moore–Penrose 广义逆矩阵，又称伪逆（pseudo inverse），记作M = A+。更进一步对复数矩阵A来说，若得出的广义逆矩阵的第 三个条件扩展为MA与AM 均为Hermit 矩阵，则这样构造的矩阵也是唯一的。 MATLAB 下B=pinv(A) 可以求出A矩阵的Moore–Penrose 广义逆矩阵B。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解77 3.1.2 矩阵的分解 矩阵的相似变换在状态空间分析中有着重要的意义。这里将介绍与矩阵变换 与分解方面的内容，如三角分解、奇异值分解等。 1. 矩阵的相似变换 假设有一个n × n的方阵A，并存在一个和它同阶的非奇异矩阵T ，则可以对 A矩阵进行如下的变换 b A = T ?1AT （3-1-5） 这种变换称为A的相似变换（similarity transform）。可以证明，变换后矩阵b A的特 征值和原矩阵A是一致的，亦即相似变换并不改变原矩阵的特征结构。 2. 矩阵的三角分解 矩阵的三角分解又称为LU 分解，其目的是将一个矩阵分解成一个下三角矩阵 L和一个上三角矩阵U 的乘积，即A = LU，其中L和U 矩阵可以分别写成 L =2664 1 l21 1 ... ... ... ln1 ln2 · · · 1 3775 , U =2664 u11 u12 · · · u1n u22 · · · u2n ... ... unn 3775 （3-1-6） MATLAB 下提供了[L,U]=lu(A) 函数，可以对给定矩阵A进行LU 分解。 如果采用数值算法，则主元素可能进行必要的换行，所以有时得出的L矩阵不是真 正的下三角矩阵，而是其基本置换形式；如果A为符号型矩阵，则可以得出真正的 三角分解。 3. 对称矩阵的Cholesky分解 如若A矩阵为对称矩阵，则仍然可以用LU 分解的方法对之进行分解，对称矩 阵LU 分解有特殊的性质，即L=UT，令DT=L为一个下三角矩阵，则可以将原来 矩阵A分解成A=DTD，其中D矩阵可以形象地理解为原A矩阵的平方根。对该 对称矩阵进行分解可以采用Cholesky 分解算法。MATLAB 提供了chol() 函数来 求取矩阵的Cholesky 分解矩阵D，该函数的调用格式可以写成[D,p]=chol(A)， 式中，返回的D为Cholesky 分解矩阵，且A=DTD；而p ? 1 为A矩阵中正定的子 矩阵的阶次，如果A为正定矩阵，则返回p = 0。 4. 矩阵的正交基 对于一类特殊的相似变换矩阵T 来说，如果它本身满足T ?1 = T ?，其中T ? 为T 的Hermit 共轭转置矩阵，则称T 为正交矩阵，并将之记为Q = T 。可见正交矩 阵Q满足下面的条件： Q?Q = I， 且QQ? = I （3-1-7） 其中I 为n × n的单位矩阵。MATLAB 中提供了Q=orth(A) 函数来求A矩阵的 正交基Q，其中Q的列数即为A矩阵的秩。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 5. 矩阵的奇异值分解 假设A矩阵为n×m矩阵，且rank(A) = r，则A矩阵可以分解为A = LΛMT， 其中L和M均为正交矩阵，Λ = diag(σ1, σ2, · · · , σn) 为对角矩阵，其对角元素 σ1, σ2, · · · , σn 满足不等式σ1 ? σ2 ? · · · ? σn ? 0。 MATLAB 提供了直接求矩阵奇异值分解的函数[L,A1,M]=svd(A)，其中， A为原始矩阵，返回的A1 为对角矩阵，而L和M均为正交变换矩阵，并满足A = LA1MT。 6. 矩阵的条件数 矩阵的奇异值大小通常决定矩阵的性态，如果矩阵的奇异值的差异特别大，则 矩阵中某个元素有一个微小的变化将严重影响到原矩阵的参数，这样的矩阵又称 为病态矩阵或坏条件矩阵，而在矩阵存在等于0 的奇异值时称为奇异矩阵。矩阵 最大奇异值σmax 和最小奇异值σmin 的比值又称为该矩阵的条件数，记作cond(A)， 即cond(A) = σmax/σmin，矩阵的条件数越大，则对元素变化越敏感。矩阵的最大和 最小奇异值还分别经常记作ˉσ(A) 和σ(A)。在MATLAB 下也提供了函数cond(A) 来求取矩阵A的条件数。 3.1.3 矩阵指数eA 和指数函数eAt 如果已知方阵A，则该矩阵的指数定义为 eA = ∞Xi=0 1 i! Ai = I + A + 1 2! A2 + 1 3! A3 + · · · + 1 m! Am + · · · （3-1-8） 矩阵指数可由MATLAB 给出的expm(A) 函数立即求出，矩阵的其他函数，如 cosA可以由funm(A,@cos) 函数求出。值得指出的是：funm() 函数采用了特征值、 特征向量的求解方式。若矩阵含有重特征根，则特征向量矩阵为奇异矩阵，这样该 函数将失效，这时应该考虑用Taylor 幂级数展开的方式进行求解[6]。一般矩阵函数 还可以考虑文献[1]中介绍的解析解方法。 例3-3 已知矩阵如下，试求该矩阵的指数与指数函数。 A = 24 ?11 ?5 5 12 5 ?6 0 1 035 解矩阵指数eA 和指数函数eAt 可以由下面语句直接求出： >> A=[-11,-5,5; 12,5,-6; 0,1,0]; expm(A) % 求数值解 A=sym(A); expm(A), syms t; expm(A*t) % 求解析解和指数函数 指数矩阵的数值解、解析解与eAt 分别为 eA ≈ 24 0.24737701 0.30723864 0.42774107 0.14460292 ?0.00080693 ?0.51328929 0.88197566 0.82052793 0.3064317135 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解79 eA =24 15e?3 ? 20e?2 + 6e?1 5e?1 ? 15e?2 + 10e?3 5e?2 ? 5e?3 24e?2 ? 18e?3 ? 6e?1 ?12e?3 ? 5e?1 + 18e?2 ?6e?2 + 6e?3 6e?1 ? 12e?2 + 6e?3 ?9e?2 + 4e?3 + 5e?1 ?2e?3 + 3e?2 35 eAt=24 15e?3t ? 20e?2t + 6e?t 5e?t ? 15e?2t + 10e?3t 5e?2t ? 5e?3t 24e?2t ? 18e?3t ? 6e?t ?12e?3t ? 5e?t + 18e?2t ?6e?2t + 6e?3t 6e?t ? 12e?2t + 6e?3t ?9e?2t + 4e?3t + 5e?t ?2e?3t + 3e?2t35 3.1.4 矩阵的任意函数计算 矩阵任意函数f(A) 的数学定义为 f(A) = I+f(0)A+ 1 2! f˙(0)A2+ 1 3! ¨ f(0)A3+· · ·+ 1 m! f(m?1)(0)Am+· · · （3-1-9） 很多矩阵函数可以由funm() 直接求解，也可以利用文献[1]给出的funmsym() 函数直接求取。下面将通过例子演示该函数的调用方法。 例3-4 已知矩阵A如下，试求复合矩阵函数ψ(A) = eAcosAt。 A = 264 ?7 2 0 ?1 1 ?4 2 1 2 ?1 ?6 ?1 ?1 ?1 0 ?4 375 解如果想求出ψ(A) = eAcosAt，则应该构造原型函数为f=exp(x*cos(x*t))，这样 就可以用下面语句直接求取矩阵函数了。 >> A=[-7,2,0,-1; 1,-4,2,1; 2,-1,-6,-1; -1,-1,0,-4]; % 输入矩阵 syms x t; A=sym(A); A1=funmsym(A,exp(x*cos(x*t)),x) % 直接运算 A2=expm(A*funm(A*t,@cos)) % 新版本MATLAB 下还可以使用这个命令 得出的结果是很冗长的，这里不给出显示的内容。 3.2 代数方程的MATLAB求解 方程求解是科学与工程领域到处都会遇到的问题。本节首先介绍各类线性方 程的解析解与数值解方法，然后计算非线性方程的设置求解方法，并介绍矩阵方 程的数值解法，还将试图得出多解矩阵方程全部的根。 3.2.1 线性方程求解问题及MATLAB实现 本节将介绍各种矩阵方程的求解方法，首先介绍矩阵逆和伪逆的求解方法，然 后介绍一般线性代数方程的求解、Lyapunov 方程与Riccati 方程求解问题。 1. 线性方程求解 前面已经介绍过矩阵的左除和右除，可以用来求解线性方程。若线性方程为 AX = B，则用X=A\B即可求出方程的解；若方程为XA = B，则用X=B/A 即可求出方程的解。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 更严格地，求解线性代数方程AX = B应该分下面几种情况考虑[1]： （1）若矩阵A为非奇异方阵，则方程的唯一解为X=inv(A)*B。 （2）若A为奇异方阵，如果A和C = A,B矩阵的秩均为m，则线性代数 方程有无穷多解，这时可以由?x=null(A) 得出齐次方程Ax = 0的基础解系，用 x0=pinv(A)*B求出原方程的一个特解，这时定义符号变量a1, a2, · · · , an?m，则 原方程的解为 x=a1*?x(:,1)+a2*?x(:,2)+· · ·+an?m*?x(:,n ? m)+x0 （3）若A和C=A,B矩阵的秩不同，则方程没有解，只能用x=pinv(A)*B 命令求出方程的最小二乘解。 另外，采用rref(C) 可以对原方程进行基本行变换，得出方程的解析解。 例3-5 试求解线性代数方程组并验证得出的结果。 264 1 2 3 4 2 2 1 1 2 4 6 8 4 4 2 2 375 X = 264 1 3 2 6 375 解由下面的语句可以求出矩阵A和判定矩阵C = [A,B] 的秩。 >> A=[1 2 3 4; 2 2 1 1; 2 4 6 8; 4 4 2 2]; B=[1;3;2;6]; C=[A B]; [rank(A), rank(C)] 通过检验秩的方法得出矩阵A和C 的秩相同，都等于2，小于矩阵的阶次4，由此可 以得出结论，原线性代数方程组有无穷多组解。如需求解原代数方程组，可以先求出化 零空间Z，并得出满足方程的一个特解x0。 >> syms a1 a2; Z=null(sym(A)); x0=sym(pinv(A))*B; x=a1*Z(:,1)+a2*Z(:,2)+x0, A*x-B 由上面结果可以写出方程的解析解如下，经检验可见误差矩阵为零。 x = α1 264 2 ?5/2 1 0 375 + α2 264 3 ?7/2 0 1 375 + 264 125/131 96/131 ?10/131 ?39/131 375 = 264 2a1 + 3a2 + 125/131 ?5a1/2 ? 7a2/2 + 96/131 a1 ? 10/131 a2 ? 39/131 375 如果采用D=rref(C) 函数，利用基本行变换得出方程的解析解，得出 D = 264 1 0 ?2 ?3 2 0 1 5/2 7/2 ?1/2 0 0 0 0 0 0 0 0 0 0 375 这样可以写出方程的解为x1 = 2x3 + 3x4 + 2，x2 = ?5x3/2 ? 7x4/2 ? 1/2，其中， x3，x4 可以取任意常数。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解81 2.AXB = C 型方程 如果已知矩阵为相容矩阵，则可以利用Kronecker 乘积技术将原方程变换为 ??BT ? Ax = c （3-2-1） 其中，?为Kronecker 乘积，后面将给出定义与求解方法，x = vec(X)，c = vec(C) 为矩阵X，C 按列展开的列向量。其中，MATLAB 命令C(:) 可以将矩阵C 按列展 开，获得c 列向量，而由C=reshape(c,n,m) 命令则可以从向量变换回矩阵C，如 果维数n、m已知；若维数未知，则可以由C1=reshape(c,size(C)) 命令还原。 矩阵A、B的Kronecker 乘积A ? B的数学定义为 A ? B =26664 a11B a12B · · · a1nB a21B a22B · · · a2nB ... ... ... ... am1B am2B · · · amnB 37775 （3-2-2） 矩阵的Kronecker 乘积可以由kron() 函数直接计算，C=kron(A,B)。 例3-6 试求解如下的矩阵方程。 24 8 1 6 3 5 7 4 9 235 X26664 0 1 0 0 1 1 0 1 2 2 1 2 0 0 2 0 0 1 1 1 1 0 0 2 1 37775 =24 0 2 0 0 2 1 2 1 0 0 2 1 1 1 035 解可以利用下面的语句直接求解给出的方程，得出方程的解析解。 >> B=[0,1,0,0,1; 1,0,1,2,2; 1,2,0,0,2; 0,0,1,1,1; 1,0,0,2,1]; C=[0,2,0,0,2; 1,2,1,0,0; 2,1,1,1,0]; A=[8,1,6; 3,5,7; 4,9,2]; x=inv(kron(sym(B'),A))*C(:) X=reshape(x,size(C)), A*X*B-C 得出的结果如下，经验证该解确实满足原始方程，误差为零。 X =24 257/360 7/15 ?29/90 ?197/360 ?29/180 ?179/180 ?8/15 23/45 119/180 23/90 ?163/360 ?8/15 31/90 223/360 31/18035 3. Lyapunov方程求解 下面的方程称为Lyapunov 方程 AX +XAT = ?C （3-2-3） 其中A、C 为给定矩阵，且C 为对称矩阵。MATLAB 下提供的X=lyap(A,C) 可 以立即求出满足Lyapunov 方程的矩阵X。该函数亦可用于不对称C 矩阵时方程 的求解。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 描述离散系统的Lyapunov 方程标准型为 AXAT ?X + Q = 0 （3-2-4） 该方程可以直接用MATLAB 现成函数dlyap() 求解，即X=dlyap(A,Q)。 4. Sylvester方程求解 Sylvester 方程实际上是Lyapunov 方程的推广，有时又称为广义Lyapunov 方 程，该方程的数学表示为 AX +XB = ?C （3-2-5） 其中A、B、C 为给定矩阵。MATLAB 下提供的X=lyap(A,B,C) 可以立即求出 满足该方程的X矩阵。 其实可以证明，利用Kronecker 乘积可以将方程改写成线性代数方程。 ??Im ? A + BT ? Inx = ?c （3-2-6） 其中，c = vec(C)，x = vec(X) 为矩阵按列展开的列向量。 文献[1,3]基于该方法给出了求取一般Lyapunov 方程和Sylvester 方程的解析 解函数lyapsym()。针对不同的方程类型，可以由下面的格式分别求解。 X=lyapsym(sym(A),C), % 连续Lyapunov 方程 X=lyapsym(sym(A),-inv(A'),Q*inv(A')), % 离散Lyapunov 方程 X=lyapsym(sym(A),B,C), % Sylvester 方程 例3-7 求解下面的Sylvester 方程。 24 8 1 6 3 5 7 4 9 2 35 X +X 24 16 4 1 9 3 1 4 2 1 35 = 24 1 2 3 4 5 6 7 8 0 35 解调用lyap() 函数可以立即得出原方程的数值解。 >> A=[8,1,6; 3,5,7; 4,9,2]; B=[16,4,1; 9,3,1; 4,2,1]; C=-[1,2,3; 4,5,6; 7,8,0]; X=lyap(A,B,C), norm(A*X+X*B+C) 可以得出该方程的数值解如下，经检验该解的误差为9.5337×10?15，精度较高。 X = 24 0.0749 0.0899 ?0.4329 0.0081 0.4814 ?0.2160 0.0196 0.1826 1.157935 如果想获得原方程的解析解，则可以使用下面的语句直接求解。 >> x=lyapsym(sym(A),B,C), norm(A*x+x*B+C) 得出方程的解如下，经检验误差为零，该解是原方程的解析解。 x = 264 1349214/18020305 648107/7208122 ?15602701/36040610 290907/36040610 3470291/7208122 ?3892997/18020305 70557/3604061 1316519/7208122 8346439/7208122 375 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解83 5. Riccati方程求解 下面的方程称为Riccati 代数方程。 ATX +XA ?XBX + C = 0 （3-2-7） 其中A、B、C 为给定矩阵，且B为非负定对称矩阵，C 为对称矩阵，则可以通过 MATLAB 的are() 函数得出Riccati 方程的解：X=are(A,B,C)，且X为对称 矩阵。离散系统的Riccati 方程可以用dare() 函数直接求解。 例3-8 试求解并检验下面给出的Riccati 方程。 24 ?2 ?1 0 1 0 ?1 ?3 ?2 ?235 X +X24 ?2 1 ?3 ?1 0 ?2 0 ?1 ?235 ?X24 2 2 ?2 ?1 5 ?2 ?1 1 235 X +24 5 ?4 4 1 0 4 1 ?1 535 = 0 解对比所述方程和式（3-2-7）给出的标准型可见 A =24 ?2 1 ?3 ?1 0 ?2 0 ?1 ?235 , B =24 2 2 ?2 ?1 5 ?2 ?1 1 235 , C =24 5 ?4 4 1 0 4 1 ?1 535 可以用下面的语句直接求解该方程，经验证得出解的误差为1.4215×10?14。 >> A=[-2,1,-3; -1,0,-2; 0,-1,-2]; B=[2,2,-2; -1 5 -2; -1 1 2]; C=[5 -4 4; 1 0 4; 1 -1 5]; X=are(A,B,C); norm(A'*X+X*A-X*B*X+C) 原方程的数值解为 X =24 0.98739 ?0.798330 0.41887 0.57741 ?0.130790 0.57755 ?0.28405 ?0.073037 0.6924135 3.2.2 一般非线性方程的求解 非线性方程如果不借助计算机也是难以求解的问题。本节将探讨一般非线性 方程的准解析解方法，还将介绍二元非线性方程的图解方法与一般非线性方程的 数值解方法。 1. 非线性方程的解析解法 MATLAB 符号运算工具箱中提供了solve() 函数可以直接求出某些方程的 解析解。如果方程没有解析解，还可以通过vpasolve() 函数得出方程的高精度数 值解。用户只需用符号表达式表示出所需求解的方程即可以直接得出方程的解析 解或高精度数值解。值得指出的是，早期版本允许使用字符串描述方程，不过在新 版本下这样的描述已经不再支持了，只能用符号表达式描述方程。该函数尤其适用 于可以转化成多项式类方程的准解析解[1]。 例3-9 试求解鸡兔同笼问题，其数学形式为下面给出的二元一次方程组。 ( x + y = 35 2x + 4y = 94 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 解求解鸡兔同笼问题有各种各样的方法，这里介绍基于vpasolve() 函数的求解方法， 需要首先将方程表示为符号表达式，然后直接求解方程即可。注意，方程中的等号需要 用== 表示；如果方程右侧为零，则可以略去==0。 >> syms x y; [x0,y0]=vpasolve(x+y==35,2*x+4*y==94) 例3-10 试求解下面给出的联立方程。 8> <> : 1 2 x2 + x + 3 2 + 2 1 y + 5 2y2 + 3 1 x3 = 0 y 2 + 3 2x + 1 x4 + 5y4 = 0 解用手工方法不可能求解原方程，所以应该考虑采用fsolve() 函数直接求解。这样由 下面的MATLAB 语句可以得出原方程全部26 个根，全部为共轭复数根。可以看出，求 解这样复杂的方程对用户而言，和求解鸡兔同笼问题一样简单。 >> syms x y [x,y]=vpasolve(x^2/2+x+3/2+2/y+5/(2*y^2)+3/x^3==0,... y/2+3/(2*x)+1/x^4+5*y^4==0); size(x) 2. 非线性方程的图解法 前面介绍过，满足隐式方程的解可以由fimplicit() 函数或ezplot() 函数直 接绘制出来。如果想求出若干个隐式方程构成的联立方程的解，则可以将这些隐式 方程用fimplicit() 函数或ezplot() 在同一坐标系下绘制出来，这样，这些曲线 的交点就是原联立方程的解，可以利用局部放大的方法把感兴趣的解从图形上读 出来。具体求解方法可以参见例2-36 中给出的演示。 3. 一般非线性方程的MATLAB数值解法 前面介绍了非线性方程组的两种解法，但这些解法均有一定的局限性。例如， solve() 函数适合于求解可以转换成多项式形式的方程解，对一般超越方程没有 较好的解决方法，而图解法适合求解一元、二元方程的解，且求解精度由于坐标轴 数据显示只能保留小数点位数有限，所以精度较低，且只能得到方程的实根。 MATLAB 提供了fsolve() 函数，利用搜索的方法求解一般非线性方程组。该 函数求解一般非线性方程组的步骤如下： （1）变换成方程的标准型。Y =F(X)=0，其中，X和F 是同维数的矩阵。 （2）用MATLAB描述方程。可以采用匿名函数或M-函数直接描述方程。 （3）选定初值求解方程。求解函数调用格式为 [x,f1,flag,details]=fsolve(fun,x0,options) 其中，fun 为步骤（2）中建立的方程组MATLAB 表示，x0 为给定的初值。变量 x为搜索出来的方程的解，f1 为该解带入原方程得出的误差值。返回的flag 变量 为标志量，如果其值大于0 表示求解成功。details 还将返回一些中间信息，如迭 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解85 代步数等。如果想修改求解的误差限，则可以设置options 模板，该值可以通过 optimset() 函数设定。 （4）解的验证。将得出的解代入方程F=fun(X)，由norm(F) 检验结果。 例3-11 试利用数值解法重新求解例2-36 中给出的超越方程。 解由于该方程是关于自变量x 和y 的，而标准型方程是针对自变量x的，所以可以考虑 引入变量x1 = x，x2 = y，这样原方程可以写成下面的标准型。 y = f(x) = "x21 e?x1x2 2/2 + e?x1/2 sin(x1x2) x22 cos(x2 + x21 ) + x21 ex1+x2 #= 0 其中，x = [x1, x2]T。这样可以由如下的匿名函数描述原始的非线性方程组。 >> f=@(x)[x(1)^2*exp(-x(1)*x(2)^2/2)+exp(-x(1)/2)*sin(x(1)*x(2)); x(2)^2*cos(x(2)+x(1)^2)+x(1)^2*exp(x(1)+x(2))]; 选定例子中得出的某随机点作为初始搜索点，则得出的解为x = 2.7800，y = ?3.3902，代入原方程可见误差达10?11 级别，可见求解精度大大增加。 >> x0=-5+10*rand(2,1); x=fsolve(f,x0); y=x(2), x=x(1) 改变初值x0，则可以得到方程其他的实根。例如，由下面的语句可能得出另一对根 x = 0，y = 1.5708，带入原方程可见精度达到10?7 级，基本满足一般要求。反复使用上 述的语句还可能得出很多其他的根。 >> x0=rand(2,1); x=fsolve(f,x0), f(x) 利用数值求解函数fsolve()，还可以人为设置求解精度等控制量，例如可以用下面 的语句直接求解方程，得到精确些的解。例如，用下面的语句重新求解原方程，则上述第 一个根的精度可以增加到10?14 级。 >> x0=[2.7795; -3.3911]; ff=optimset; ff.TolX=1e-20; ff.TolFun=1e-20; x=fsolve(f,x0,ff), f(x) 如果选择复数初值，如x0 = [5 + 4j,?1 + 2j]，则可以用下面语句求解方程，得出 x = [3.5470 + 0.0480j,?3.5422 + 0.0479j]，误差也为10?14 级。 >> x0=[5+4i; -1+2i]; x=fsolve(f,x0,ff), f(x) 3.2.3 非线性矩阵方程的MATLAB求解 前面介绍的fsolve() 函数可以直接用于求解非线性矩阵方程的一个根。如果 给定初值，则可以通过这个初值搜索出其他的根。若给出多个初值，则可能求出其 他的根。可以编写一个求解函数，一次性得出方程的多个根，该函数是文献[4]给出 函数的改进形式，方程运行时间越长得出的结果可能越精确。 function more_sols(f,X0,varargin) [A,tol,tlim,ff]=default_vals({1000,eps,30,optimset},varargin{:}); if length(A)==1, a=-0.5*A; b=0.5*A; else, a=A(1); b=A(2); end ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. ar=real(a); br=real(b); ai=imag(a); bi=imag(b); ff.Display='off'; ff.TolX=tol; ff.TolFun=1e-20; [n,m,i]=size(X0); X=X0; tic if i==0, X0=zeros(n,m); % 判定零矩阵是不是方程的孤立解 if norm(f(X0))1e-5, x0=x0+(ai+(bi-ai)*rand(n,m))*1i; end [x,aa,key]=fsolve(f,x0,ff); t=toc; if t>tlim, break; end if key>0, N=size(X,3); % 读出已记录根的个数，若该根已记录，则放弃 for j=1:N, if norm(X(:,:,j)-x)<1e-4; key=0; break; end, end if key==0 % 如果找到的解比存储的更精确，则替换 if norm(f(x))0, X(:,:,i+1)=x; % 记录找到的根 if norm(imag(x))>1e-5 && norm(f(conj(x)))<1e-5 i=i+1; X(:,:,i+1)=conj(x); % 若找到复根，则测试其共轭复数 end assignin('base','X',X); i=i+1, tic % 更新信息 end, assignin('base','X',X); end, end 其子函数default_val() 的清单为 function varargout=default_vals(vals,varargin) if nargout~=length(vals), error('number of arguments mismatch'); else, nn=length(varargin)+1; varargout=varargin; for i=nn:nargout, varargout{i}=vals{i}; end, end, end 该函数调用格式为more_sols(f,X0,A,?,tlim)，其中，f 为原函数的MATLAB 表示，可以为匿名函数、MATLAB 函数等，其他的参数可以采用默认的值，一 般无经验的用户无须给出。X0 是一个三维数组，表示以前已经得到的方程根，A 为随机数初值范围，表示初值在(?A/2,A/2) 范围内取均匀分布随机数，默认值为 1000。? 为求解的默认精度要求，默认值为5 倍的eps。tlim 为允许的等待时间，默认 值为30，表示半分钟，即半分钟内如果没有找到新解，则停止搜索。 观察图2-13（a）中的曲线交点，可以看出，x = 0，y = 0 这个点不是由两条曲线 演化而来，所以这类解称为方程的孤立根，是不同通过搜索方法得出的。在上面给 出的求解函数中，首先尝试原点是不是孤立解，如果是则记录下来。另外，如果搜索 到了一个复数根，则尝试一下其共轭复数是不是方程的根，如果是也记录下来，这 样可以通过程序的使用效率。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解87 由于使用了死循环while(1)，只能由用户给出的中断命令Ctrl+C 键停止运 行，或等待tlim 后没有新解后自动终止，这样该函数不能返回任何变元。为解决这 样的问题，使用assignin() 函数将得到的解X和求解方程次数M 写入MATLAB 的工作空间。其中，X为三维数组，X(:,:,i) 对应于第i 个根。已找到根的个数可 以由size(X,3) 读出。 例3-12 试求出例3-8 中给出的Riccati 方程全部的根。 解由于该方程是关于X的二次型方程，从常理看该方程可能存在其他的根，但are() 函 数只能求出一个根。其他的根可以使用搜索的方法求出。现在可以试用more_sols() 函 数来求解Riccati 方程其他可能的根。首先，将矩阵方程用下面匿名函数直接描述出来。 >> A=[-2,1,-3; -1,0,-2; 0,-1,-2]; B=[2,2,-2; -1 5 -2; -1 1 2]; C=[5 -4 4; 1 0 4; 1 -1 5]; f=@(X)A'*X+X*A-X*B*X+C; more_sols(f,zeros(3,3,0)) 可见，上述函数的定义格式和Riccati 方程的数学描述一样简洁。定义了方程函数， 则直接调用more_sols() 函数即可得出方程的解。该方程有8 个根，所以显示i = 8 以 后就可以用Ctrl+C 键结束程序运行，或等待程序自动停止。这时，工作空间中的X三维 数组将返回方程的所有8 个根，前面are() 函数求出的解矩阵只是其中之一。 X1 =24 0.8878 ?0.9608 ?0.2446 0.1071 ?0.8984 ?2.5562 ?0.0185 0.3604 2.4619 35 , X2 =24 ?0.1538 0.1086 0.4622 2.0277 ?1.7436 1.3474 1.9003 ?1.7512 0.5057 35 X3 =24 1.2212 ?0.4165 1.9775 0.3577 ?0.4893 ?0.8863 ?0.7414 ?0.8197 ?2.3559 35 , X4 =24 ?2.1032 1.2977 ?1.9697 ?0.2466 ?0.3563 ?1.4899 ?2.1493 0.7189 ?4.5464 35 X5 =24 0.9873 ?0.7983 0.4188 0.5774 ?0.1307 0.5775 ?0.2840 ?0.0730 0.6924 35 , X6 =24 0.6664 ?1.3222 ?1.7200 0.3120 ?0.5640 ?1.1910 ?1.2272 ?1.6129 ?5.5939 35 X7 =24 ?0.7618 1.3312 ?0.8400 1.3182 ?0.3173 ?0.1718 0.6371 0.7884 ?2.1996 35 , X8 =24 23.947 ?20.667 2.4528 30.146 ?25.983 3.6699 51.967 ?44.911 4.6409 35 如果将求解范围用复数表示，则还可以得出方程的复数根。可见，通过搜索可以得 出方程全部20 个根。 >> more_sols(f,X,1000+1000i) 利用这里给出的more_sols() 函数可以求解由其他方法很难求解的矩阵方程， 例如下面Riccati 方程扩展形式得出的特殊方程也可以直接求解。 AX +XD ?XBX + C = 0 （3-2-8） AX +XD ?XBXT + C = 0 （3-2-9） ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 例3-13 试求解式（3-2-9）中给出的矩阵方程，其中 A = 24 2 1 9 9 7 9 6 5 3 35 , B = 24 0 3 6 8 2 0 8 2 8 35 , C = 24 7 0 3 5 6 4 1 4 4 35 , D = 24 3 9 5 1 2 9 3 3 0 35 解可以通过下面的语句直接求出该方程的16 个实根，从略。 >> A=[2,1,9; 9,7,9; 6,5,3]; B=[0,3,6; 8,2,0; 8,2,8]; C=[7,0,3; 5,6,4; 1,4,4]; D=[3,9,5; 1,2,9; 3,3,0]; f=@(X)A*X+X*D-X*B*X.'+C; more_sols(f,zeros(3,3,0)) 例3-14 例2-36 中的方程也可以看成是一种特殊的矩阵方程，试求解该方程。 解仿照例3-11 中介绍的方法，由匿名函数描述该方程，则可以给出如下的命令，直接求 解该方程。 >> f=@(x)[x(1)^2*exp(-x(1)*x(2)^2/2)+exp(-x(1)/2)*sin(x(1)*x(2)); x(2)^2*cos(x(2)+x(1)^2)+x(1)^2*exp(x(1)+x(2))]; more_sols(f,zeros(2,1,0),4*pi) 经过一段时间的等待，或用Ctrl+C 键强制停止搜索过程，用下面的命令将得出的解 叠印在图解法得出的图上，如图3-1 所示，搜索到的解用圈表示。可见，绝大部分的交点 均已找到。将所有的解代入原方程，则可以得出最大误差e = 5.9706×10?13，可见得出 的解精度远远高于图解法。如果程序运行一段时间后自然停止，则可能找到该方程在感 兴趣区域的全部41 个根。 --6 --4 --2 0 2 4 6 --6 --4 --2 0 2 4 6 图3-1 联立方程图解法及搜索到的解 >> syms x y; f1=x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y); f2=y^2*cos(y+x^2)+x^2*exp(x+y); fimplicit([f1,f2],[-2*pi,2*pi],'MeshDensity',800), hold on x0=X(1,1,:); x0=x0(:); y0=X(2,1,:); y0=y0(:); plot(x0,y0,'o') F=[]; for i=1:length(x0), F=[F,f([x0(i),y0(i)])]; end, norm(F) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解89 3.3 常微分方程问题的MATLAB求解 微分方程问题是动态系统仿真的核心，由强大的MATLAB 语言可以对一阶微 分方程组求取数值解，其他类型的微分方程可以通过合适的算法变换成可解的一 阶微分方程组进行求解，这里将介绍微分方程的求解方法。 3.3.1 一阶常微分方程组的数值解法 假设一阶常微分方程组由下式给出 x˙ i(t) = fi(t, x(t)), i = 1, 2, · · · , n （3-3-1） 其中，x(t) 为状态变量xi(t) 构成的向量，即x(t) = [x1(t), x2(t), · · · , xn(t)]T，常称 为系统的状态向量，n称为系统的阶次，fi(·) 为任意非线性函数，t 为时间变量。一般 来说，非线性的微分方程是没有解析解的，只能采用数值方法，在初值x(0) 下求解 常微分方程组。 求解常微分方程组的数值方法是多种多样的，如常用的Euler 法、Runge–Kutta 算法、Adams 线性多步法、Gear 算法等。为解决刚性（stiff）问题又有若干专用的刚 性问题求解算法；另外，如需要求解隐式常微分方程组和含有代数约束的微分代数 方程组时，则需要对方程进行相应的变换，方能进行求解。微分方程的数值解法详 参文献[5]，本节只给出简单微分方程的直接数值求解方法。 MATLAB 中给出了若干求解一阶常微分方程组的函数，如ode23()（二阶三 级Runge–Kutta 算法）、ode45()（四阶五级Runge–Kutta 算法）、ode15s()（变阶 次刚性方程求解算法）等，其调用格式都是一致的： [t, x]=ode45(Fun,tspan,x0,options, 附加参数) 其中，t 为自变量构成的向量，一般采用变步长算法，返回的x是一个矩阵，其列 数为n，即微分方程的阶次，行数等于t 的行数，每一行对应于相应时间点处的状 态变量向量的转置。Fun 为用MATLAB 编写的固定格式的M-函数或匿名函数，描 述一阶微分方程组，tspan 为数值解时的初始和终止时间等信息，x0 为初始状态变 量，options 为求解微分方程的一些控制参数，还可以将一些附加参数在求解函数 和方程描述函数之间传递。下面将通过例子介绍微分方程求解过程。 例3-15 试求解下面给出的著名的R?ssler 微分方程组，选定a = b = 0.2, c = 5.7，且 x(0)=y(0)=z(0)=0。8> <> : x˙ (t) = ?y(t) ? z(t) y˙(t) = x(t) + ay(t) z˙(t) = b + [x(t) ? c]z(t) 解由于该方程是非线性微分方程，一般没有解析解，只能通过数值解的方法来研究该 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 方程。引入新状态变量x1(t)=x(t)，x2(t)=y(t)，x3(t)=z(t)，则可将原微分方程改写成 8> <> : x˙ 1(t) = ?x2(t) ? x3(t) x˙ 2(t) = x1(t) + ax2(t) x˙ 3(t) = b + [x1(t) ? c]x3(t) 其矩阵形式为x˙ (t) = 264 ?x2(t) ? x3(t) x1(t) + ax2(t) b + [x1(t) ? c]x3(t) 375 若想求解这个微分方程，则需要用户自己编写一个MATLAB 函数描述它。 function dx=rossler(t,x) % 虽然不显含时间，还应该写出占位 dx=[-x(2)-x(3); % 对应方程第一行，直接将参数代入 x(1)+0.2*x(2); 0.2+(x(1)-5.7)*x(3)]; % 其余两行 对比此函数和给出的数学方程，应该能看出编写这样的函数还是很直观的。只要能 得出一阶微分方程组，则可以立即编写出MATLAB 函数来描述它。编写了该函数，就可 以将其存成rossler.m 文件。除了用MATLAB 函数描述微分方程之外，对简单问题还可 以用匿名函数的方式描述原始的微分方程。 >> f=@(t,x)[-x(2)-x(3); x(1)+0.2*x(2); 0.2+(x(1)-5.7)*x(3)]; 这样就可以用下面的MATLAB 语句求出微分方程的数值解。 >> x0=[0; 0; 0]; [t,y]=ode45(@rossler,[0,100],x0); % 解方程 % 或采用[t,y]=ode45(f,[0,100],x0) 命令求解方程 plot(t,y) % 绘制各个状态变量的时间响应 figure; plot3(y(:,1),y(:,2),y(:,3)), grid % 绘制相空间轨迹 上面的命令将直接得出该微分方程在t ∈ [0, 100] 内的数值解，该数值解可以用图 形更直观地表示出来，若绘制各个状态变量和时间之间的关系曲线，则得出如图3-2（a） 所示的时域响应曲线，该方程研究的另外一种实用的表示是三维曲线绘制，用三维图形 可以绘制出相空间曲线，如图3-2（b）所示。如果用comet3() 函数取代plot3()，则可以 看出轨迹走行的动画效果。 0 20 40 60 80 100 --10 0 10 20 30 （a）状态变量的时间曲线 0 10 10 10 0 5 20 0 --10 --10 --5 （b）系统响应的相空间表示 图3-2 R?ssler 方程的数值解表示 现在演示附加参数的使用方法，假设a、b、c 这三个参数需要用外部命令给出，则可 以按下面的格式写出一个新的M-函数来描述微分方程组。 function dx=rossler1(t,x,a,b,c) % 加入附加参数 dx=[-x(2)-x(3); x(1)+a*x(2); b+(x(1)-c)*x(3)]; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解91 这样就可以用下面的语句直接求解该方程并绘制曲线了。 >> a=0.2; b=0.2; c=5.7; % 从函数外部定义这三个变量，无须修改函数本身 [t,y]=ode45(@rossler1,[0,100],x0,[],a,b,c); % 用附加参数解方程 这样编写M-函数有很多好处，例如若想改变β 等参数，没有必要修改M-函数，只需 在求解该方程时将参数代入即可，这样会很方便。假设现在想研究a = 2 时的微分方程 数值解，则可以给出下面的命令直接求解。 >> a=2; [t,y]=ode45(@rossler1,[0,100],x0,[],a,b,c); 事实上，对简单的微分方程而言，用匿名函数可以避开附加变量的使用，直接描述 带有参数的微分方程。例如，下面的语句同样可以求解前面的微分方程。值得指出的是， 若想修改方程的参数，应该重新定义匿名函数。 >> a=0.2; b=0.2; c=5.7; f=@(t,x)[-x(2)-x(3); x(1)+a*x(2); b+(x(1)-c)*x(3)]; [t,y]=ode45(f,[0,100],x0); % 利用匿名函数可以直接求解 在许多领域中，经常遇到一类特殊的常微分方程，其中一些解变化缓慢，另一 些变化快，且相差较悬殊，用ode45() 函数长时间的不出结果。这类方程常常称为 刚性方程，又称为Stiff 方程。刚性问题一般不适合由ode45() 这类函数求解，而应 该采用MATLAB 求解函数ode15s()，其调用格式与ode45() 一致。 3.3.2 常微分方程的转换 MATLAB 下提供的微分方程数值解函数只能处理以一阶微分方程组形式给 出的微分方程，所以在求解之前需要先将给定的微分方程变换成一阶微分方程组， 而微分方程组的变换中需要选择一组状态变量，由于状态变量的选择是可以比较 任意的，所以一阶显式微分方程组的变换也不是唯一的。这里介绍微分方程组变换 的一般方法。 首先考虑单个高阶微分方程的处理方法，假设微分方程可以写成 f(t, y, y˙, y¨, · · · , y(n)) = 0 （3-3-2） 比较简单的状态变量选择方法是令x1 = y, x2 = y˙, · · · , xn = y(n?1)，这样显 然有x˙ 1 = x2, x˙ 2 = x3, · · · , x˙n?1 = xn，另外，求解式（3-3-2），得出y(n) 的显式表达 式，y(n) = f?(t, y, y˙, · · · , y(n?1))，这时就可以写出该微分方程对应的一阶微分方程 组为( x˙ i = xi+1, i = 1, 2, · · · , n ? 1 x˙n = f?(t, x1, x2, · · · , xn) （3-3-3） 这样，原微分方程就可以用MATLAB 提供的常微分方程求解ode45()、ode15s()等 函数直接求解了。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 再考虑高阶微分方程组的变换方法，假设已知高阶微分方程组为 ( f(t, x(t), x˙ (t), · · · , x(m?1)(t), x(m)(t), y(t), · · · , y(n?1)(t), y(n)(t)) = 0 g(t, x(t), x˙ (t), · · · , x(m?1)(t), x(m)(t), y(t), · · · , y(n?1)(t), y(n)(t)) = 0 （3-3-4） 则仍旧可以选择状态变量x1(t) = x(t)，x2(t) = x˙ (t)，· · · ，xm(t) = x(m?1)(t)， xm+1(t)=y(t)，xm+2(t)=y˙(t)，· · · ，xm+n(t)=y(n?1)(t)，并将其代入式（3-3-4），则 ( f(t, x1(t), x2(t), · · · , xm(t), x˙m(t), xm+1(t), · · · , xm+n(t), x˙m+n(t)) = 0 g(t, x1(t), x2(t), · · · , xm(t), x˙m(t), xm+1(t), · · · , xm+n(t), x˙m+n(t)) = 0 （3-3-5） 求解该方程则可以得出x˙m(t) 与x˙m+n(t)，从而得出所需的一阶微分方程组，最终 使用MATLAB 中提供的函数求解这些高阶微分方程组。 例3-16 考虑著名的Van der Pol 方程y¨(t)+[y2(t)?1]y˙(t)+y(t)=0，已知y(0)=?0.2, y˙(0) = ?0.7，试用数值方法求出的Van der Pol 方程的解。 解由给出的方程可知，因为它不是显式一阶微分方程组，所以不能直接求解，而必须先 进行转换，再进行求解。选择状态变量x1(t) = y(t)，x2(t) = y˙(t)，则原方程可以变换成 ( x˙ 1(t)=x2(t) ˙ x2(t)=?[x21 (t)?1]x2(t)?x1(t), 其矩阵形式为x˙ =" x2(t) ?[x21 (t)?1]x2(t)?x1(t) # 这样可以写出如下描述此方程的匿名函数。 >> f=@(t,x)[x(2); -(x(1)^2-1)*x(2)-x(1)]; 由选定的状态变量可知，其初值可以描述成x0=[?0.2,?0.7]T，所以该方程最终可 以由下面的语句直接求解并绘图，得出的时间响应曲线如图3-3（a）所示，相平面曲线如 图3-3（b）所示。 >> x0=[-0.2; -0.7]; tn=20; [t,x]=ode45(f,[0,tn],x0); plot(t,x) % 显示两个状态变量的时间曲线 figure; plot(x(:,1),x(:,2)) % 相平面曲线 0 5 10 15 20 --3 --2 --1 0 1 2 3 （a）状态变量的时间曲线 --3 --2 --1 0 1 2 3 --3 --2 --1 0 1 2 3 （b）系统响应的相平面表示 图3-3 Van der Pol 方程的数值解表示 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解93 3.3.3 微分方程数值解的验证 前面介绍了微分方程的求解方法，求解步骤总结如下： （1）转换成标准型。由于现有的求解函数只能求解一阶显式微分方程组，所以 需要首先将原始的微分方程手工变换成标准形式。 （2）描述标准型。用MATLAB 函数或者匿名函数描述原始微分方程。对简单 问题而言，用匿名函数是最好的选择，但如果原始微分方程较烦琐，则只能用M-函 数描述。 （3）调用求解函数。调用求解函数ode45() 直接求解，对特殊的方程需要采用 刚性微分方程求解函数，如ode15s()。 （4）解的验证。MATLAB 采用变步长算法求解微分方程，其关键的监测指标是 容许相对误差限RelTol 的设置。默认的RelTol 控制量为10?3，相当于千分之一的 误差，其值过大，所以应该采用较小的值。如果两次选择不同的RelTol 值解的结果 一致，则可以认为得出的解是正确的，否则应该试更小的控制量。另外，选择不同的 求解算法也可以验证解的正确性。如果想追求双精度数据结构下最精确的结果，可 以将RelTol 的值设置成3×10?14。这些控制选项可以用下面的语句设置。 options=odeset; options.RelTol=1e-7; 例3-17 已知Apollo 卫星的运动轨迹(x, y) 满足下面的方程[7]。 8> ><>>: ¨x(t) = 2y˙(t) + x(t) ? μ?(x(t) + μ) r3 1(t) ? μ(x(t) ? μ?) r3 2(t) y¨(t) = ?2x˙ (t) + y(t) ? μ?y(t) r3 1(t) ? μy(t) r3 2(t) 其中，μ = 1/82.45，μ? = 1 ? μ，且 r1(t) = p(x(t) + μ)2 + y2(t), r2(t) = p(x(t) ? μ?)2 + y2(t) 若已知初值为x(0) = 1.2，x˙ (0) = 0，y(0) = 0，y˙(0) = ?1.04935751，试求解该方程。 解选择一组状态变量x1(t) = x(t)，x2(t) = x˙ (t)，x3(t) = y(t)，x4(t) = y˙(t)，这样就可 以得出标准型为 8> ><>>: x˙ 1(t) = x2(t) x˙ 2(t) = 2x4(t) + x1(t) ? μ?(x1(t) + μ)/r3 1(t) ? μ(x1(t) ? μ?)/r3 2(t) x˙ 3(t) = x4(t) x˙ 4(t) = ?2x2(t) + x3(t) ? μ?x3(t)/r3 1(t) ? μx3(t)/r3 2(t) 式中，r1(t) = q(x1(t) + μ)2 + x23 (t), r2(t) = q(x1(t) ? μ?)2 + x23 (t)，且μ = 1/82.45， μ? = 1 ? μ。 有了数学模型描述，则可以立即写出其相应的MATLAB 函数如下： function dx=apolloeq(t,x) mu=1/82.45; mu1=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. r2=sqrt((x(1)-mu1)^2+x(3)^2); % 中间变量赋值，不宜采用匿名函数 dx=[x(2); 2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3; x(4); -2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3]; 可以看出，由于在描述微分方程的语句中含有中间变量（如r1(t)、r2(t)）的计算，所 以不宜采用匿名函数描述方程，M-函数则是唯一可行的描述方法。 调用ode45() 函数可以求出该方程的数值解，得出的轨迹如图3-4（a）所示。 >> x0=[1.2; 0; 0; -1.04935751]; [t,y]=ode45(@apolloeq,[0,20],x0); plot(y(:,1),y(:,3)) 得出方程的数值解后，需要如下的语句对其检验，得出的新解如图3-4（b）所示，可见，这 样得出的解和前面的解不同。再进一步减小RelTol 值得到的解没有太大的变化，所以 这样得出的解是正确的。 >> options=odeset; options.RelTol=1e-7; [t,y]=ode45(@apolloeq,[0,20],x0,options); plot(y(:,1),y(:,3)) --1.5 --1 --0.5 0 0.5 1 1.5 --1 --0.5 0 0.5 1 （a）默认的解 --1.5 --1 --0.5 0 0.5 1 1.5 --1 --0.5 0 0.5 1 （b）正确的解 图3-4 Apollo 卫星的轨迹曲线 3.3.4 线性常微分方程的解析求解 由微分方程理论可知，常系数线性微分方程是存在解析解的，变系数的线性微 分方程的可解性取决于其特征方程的可解性，一般是不可解析求解的，非线性的微 分方程是不存在解析解的，在MATLAB 语言中提供了dsolve() 函数，可以用于线 性常系数微分方程的解析解求解。求解微分方程时，首先应该用syms 命令声明符 号变量，以区别于MATLAB 语言的常规数值变量，然后就可以用dsolve(表达式) 命令直接求解了。类似于前面介绍的solve() 和vpasolve() 函数，微分方程应该 由符号表达式描述。下面通过例子来演示该函数的使用方法。 例3-18 试求解下面的常系数线性微分方程。 d4y(t) dt4 + 11 d3y(t) dt3 + 41 d2y(t) dt2 + 61 dy(t) dt + 30y(t) = e?6t cos 5t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解95 解可以采用下面的MATLAB 语句求解该微分方程。其中，引入了若干个中间变量，此 外，微分方程的符号表达式中，等号应该由双等号表示。 >> syms t y(t) % 声明符号变量 y1=diff(y); y2=diff(y1); y3=diff(y3); Y=dsolve(diff(y)+11*y3+41*y2+61*y1+30*y==exp(-6*t)*cos(5*t)); pretty(simplify(Y)) % 以更好看的形式显示解析解结果 上面的语句得出结果的可读性不是很好，这里采用LATEX 变换后的结果为 y(t) = ? 79e?6t 181220 cos 5t + 109e?6t 181220 sin 5t + C1e?t + C2e?2t + C3e?3t + C4e?5t 其中，Ci 为待定系数，应该由方程的初值或边值等求出，dsolve()函数可以直接求出 带有初值或边值的微分方程解。例如已知方程的初值边值条件为y(0) = 1，y˙(0) = 1， ¨y(0) = 0，y(3)(0) = 0，则可以由下面的语句求出方程的解。 >> Y=dsolve(diff(y)+11*y3+41*y2+61*y1+30*y==exp(-6*t)*cos(5*t),... y(0)==1,y1(0)==1,y2(0)==0,y3(0)==0); 则可以得出该微分方程的解析解为 y(t)=? 79e?6t 181220 cos 5t+ 109e?6t 181220 sin 5t+ 611 80 e?t? 1562 123 e?2t+ 921 136 e?3t? 443 624 e?5t 3.4 最优化问题的MATLAB求解 最优化方法在系统仿真与控制系统计算机辅助设计中占有很重要的地位，求 解最优化问题的数值算法有很多，MATLAB 中提供了各种各样的最优化问题求解 函数，可以求解无约束最优化问题、有约束最优化问题及线性规划、二次型规划问 题等，还实现了基于最小二乘算法的曲线拟合方法。 3.4.1 无约束最优化问题求解 无约束最优化问题的一般描述为 min x f(x) （3-4-1） 其中，x = [x1, x2, · · · , xn]T，该数学表示的含义亦即求取一个x 向量，使得标量最 优化目标函数f(x) 的值为最小，故这样的问题又称为最小化问题。其实，最小化是 最优化问题的通用描述，它不失普遍性。如果要求解最大化问题，那么只需给目标 函数f(x) 乘以?1 就能立即将原始问题转换成最小化问题。 MATLAB 提供了基于单纯形算法[8] 求解无约束最优化的fminsearch() 函数， 该函数的调用格式为：[x,fopt,key,c]=fminsearch(Fun,x0,options)，其中， Fun 为要求解问题的数学描述，它可以是一个MATLAB 函数，也可以是一个函 数句柄，x0 为自变量的起始搜索点，需要用户自己去选择，options 为最优化工具 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 箱的选项设定；x为返回的解；而fopt 是目标函数在x点处的值。返回的key 表示函 数返回的条件，1 表示已经求解出方程的解，而0 表示未搜索到方程的解。返回的c 为解的附加信息，该变元为一个结构体变量，其iterations 成员变量表示迭代的 次数，而其中的成员funcCount 是目标函数的调用次数。MATLAB 的最优化工具 箱中提供的fminunc() 函数与fminsearch() 功能和调用格式均很相似，有时求解 无约束最优化问题可以选择该函数。在第7、8 章中将介绍基于数值最优化算法的 最优控制器设计方法。 另外，如果决策变量x需要满足xm ? x ? xM 前提条件，则可以采用后面将 介绍的有约束最优化求解方法，或采用John D’Errico 编写的fminsearchbnd() 函 数直接求解[9]，第11 章将直接使用该函数求解最优控制器设计问题。 3.4.2 有约束最优化问题求解 有约束非线性最优化问题的一般描述为 min x s.t. 8> >>><>>>>: Ax?B Aeqx=Beq xm?x?xM C(x)?0 Ceq(x)=0 f(x) （3-4-2） 其中，x = [x1, x2, · · · , xn]T。在约束条件中直接给出了线性等式约束Aeqx=Beq， 线性不等式约束Ax ? B，一般非线性等式约束Ceq(x) = 0，非线性不等式约束 C(x)?0和决策变量的上下界约束xm ? x ? xM。注意，这里的不等式约束全部 是?不等式，若原问题关系为?，则可以将不等式两端同时乘以?1 就能将其转换 成?不等式。 该数学表示的含义亦即求取一组x向量，使得函数f(x) 在满足全部约束条 件的基础上最小化。而满足所有约束的问题称为可行问题（feasible problem）。 MATLAB 最优化工具箱中提供了一个fmincon() 函数，专门用于求解各种约 束下的最优化问题。该函数的调用格式为： [x,fopt,key,c]=fmincon(Fun,x0,A,B,Aeq,Beq,xm,xM,CFun,OPT) 其中，Fun 为给目标函数写的M-函数，x0 为初始搜索点。各个矩阵约束如果不存在， 则应该用空矩阵来占位。CFun 为给非线性约束函数写的M-函数，OPT 为控制选项。 最优化运算完成后，结果将在变元x中返回，最优化的目标函数将在fopt 变元中返 回，选项有时是很重要的。返回变元key 若不是正数，则说明这时因故未发现原问 题的解，可以考虑改变初值，或修改控制参数OPT，再进行寻优，以得出期望的最 优值。另外，如果发现最优化问题不是可行问题，则在求解结束后将给出提示：“No feasible solution found（未找到可行解）”。Fun 变量还可以用结构体的形式描述原 问题。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解97 例3-19 试求解下面的非线性最优化问题。 min x s.t.8>><>>: x1+x2?0 ?x1x2+x1+x2?1.5 x1x2??10 ?10?x1,x2?10 ex1(4x21 + 2x22 + 4x1x2 + 2x2 + 1) 解若想求解本最优化问题，应该用下面的语句先描述出目标函数和约束函数。 function y=c3exmobj(x) y=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); function [c,ce]=c3exmcon(x) ce=[]; c=[x(1)+x(2); x(1)*x(2)-x(1)-x(2)+1.5; -10-x(1)*x(2)]; 注意，约束函数应该返回两个变元，不等式约束c 和等式约束ce，其中第2、3 条约束 都应该先乘以?1 变换成?不等式。另外，由约束条件可知，第1 条约束实际上是线性不 等式约束，所以还可以用定义A、B矩阵的形式来描述，但这样需要从M-函数中先剔除 第1 条约束。调用非线性最优化问题求解函数可以得出如下结果。 >> A=[]; B=[]; Aeq=[]; Beq=[]; xm=[-10; -10]; xM=[10; 10]; x0=[5;5]; ff=optimset; ff.TolX=1e-10; ff.TolFun=1e-20; x=fmincon(@c3exmobj,x0,A,B,Aeq,Beq,xm,xM,@c3exmcon,ff) 该语句给出的显示为“Local minimum possible（可能是局部最小值）”。这样得出的“最 优解”为xT = [0.4195, 0.4195]，可以考虑用得出的“最优解”作为初值再进一步求解， 如此可以利用循环结构得出原问题的真正最优解为xT = [?9.5474, 1.0474]，循环次数 为i = 5。 >> i=1; x=x0; while (1) [x,a,b]=fmincon(@c3exmobj,x,A,B,Aeq,Beq,xm,xM,@c3exmcon,ff); if b>0, break; end, i=i+1; % 如果求解成功则结束循环 end 有约束最优化还有几种特殊的形式，如线性规划问题、二次型规划问题，可以 使用最优化工具箱中的linprog() 和quadprog() 函数直接求解。此外，整数规划、 0-1 规划等问题可以下载专门的工具求解[1]。 3.4.3 全局最优解的尝试 与线性规划等凸问题不同，常规的非线性规划问题可能含有局部最优解。如 果初始搜索点选择不当，则很容易陷入局部最优解。仿照前面介绍的多解代数 方程求解方法，很自然地考虑在fminsearch() 或fmincon() 函数外再加一层循 环，每步循环在允许的区域内随机生成初始搜索点，则可以构造出新的求解函数 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. fminunc_global()和fmincon_global()，与底层函数相比，这两个函数更可能获 得原始问题的全局最优解。这两个函数的调用格式为 [x,f0]=fminunc_global(fun,a,b,n,N) [x,f0]=fmincon_global(fun,a,b,n,N, 其他参数) 其中，fun 可以是结构体变量，也可以是目标函数的函数句柄；a 与b 为决策变量所 在的区间，如果xm 与xM 为有限的向量，还可以直接将a 和b 设置成xm 和xM 向量； n为决策变量的个数；N 为每次寻优中，底层函数fminsearch() 或fmincon() 的 调用次数，通常情况下N = 5 ～ 10 就足够。 文献[4]对这两个函数与MATLAB 全局优化工具箱提供遗传算法、粒子群算 法等智能优化算法进行了对比研究。对测试的例子而言，这两个函数在全局性与求 解耗时方面明显优于所比较的智能优化算法。这里只通过例子演示求解全局最优 解函数的使用方法与优势。 例3-20 试找出下面非线性规划问题的全局最优解。 min q,w,k s.t. 8> >><>>>: q3+9.625q1w+16q2w+16w2+12?4q1?q2?78w=0 16q1w+44?19q1?8q2?q3?24w=0 2.25?0.25k?q1?2.25+0.25k 1.5?0.5k?q2?1.5+0.5k 1.5?1.5k?q3?1.5+1.5k k 解从给出的最优化问题看，这里要求解的决策变量为q、w和k，而标准最优化方法只 能求解向量型决策变量，所以应该作变量替换，把需要求解的决策变量由决策变量向量 表示出来。对本例来说，可以引入x1 = q1，x2 = q2，x3 = q3，x4 = w，x5 = k，另外，需要 将一些不等式进一步处理，可以将原始问题手工改写成 min x s.t. 8> >>>>><>>>>>>: x3+9.625x1x4+16x2x4+16x2 4+12?4x1?x2?78x4=0 16x1x4+44?19x1?8x2?x3?24x4=0 ?0.25x5?x1??2.25 x1?0.25x5?2.25 ?0.5x5?x2??1.5 x2?0.5x5?1.5 ?1.5x5?x3??1.5 x3?1.5x5?1.5 x5 从手工变换后的结果看，原始问题有两个非线性等式约束，没有不等式约束，所以 可以由下面语句描述原问题的非线性约束条件。 function [c,ce]=c3mnls(x) c=[]; % 非线性约束条件，其中，不等式约束为空矩阵 ce=[x(3)+9.625*x(1)*x(4)+16*x(2)*x(4)+16*x(4)^2+12... -4*x(1)-x(2)-78*x(4); 16*x(1)*x(4)+44-19*x(1)-8*x(2)-x(3)-24*x(4)]; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解99 原模型的线性约束可以写成线性不等式的矩阵形式Ax ? b，其中 A = 2666664 ?1 0 0 0 ?0.25 1 0 0 0 ?0.25 0 ?1 0 0 ?0.5 0 1 0 0 ?0.5 0 0 ?1 0 ?1.5 0 0 1 0 ?1.5 3777775 , b = 2666664 ?2.25 2.25 ?1.5 1.5 ?1.5 1.5 3777775 该问题没有线性等式约束，也没有决策变量的下界与上界约束，所以可以将这些 约束条件用空矩阵表示，或直接采用结构体描述最优化问题，可以不用考虑这些约束 的设置。为方便起见这里采用结构体形式描述原始问题。可以看出，这里的目标函数 值为x5。采用fmincon() 函数求解，很容易得出局部最优解x5 = 1.1448。现在试用 全局最优化函数求解该问题。一般情况下，每次都能得出该问题的全局最优解为x = [2.4544, 1.9088, 2.7263, 1.3510, 0.8175]T，目标函数为0.8175，耗时0.342 s。 >> clear P; P.objective=@(x)x(5); P.nonlcon=@c3mnls; P.solver='fmincon'; P.Aineq=[-1,0,0,0,-0.25; 1,0,0,0,-0.25; 0,-1,0,0,-0.5; 0,1,0,0,-0.5; 0,0,-1,0,-1.5; 0,0,1,0,-1.5]; % 结构体描述 P.Bineq=[-2.25; 2.25; -1.5; 1.5; -1.5; 1.5]; P.options=optimset; P.x0=rand(5,1); tic, [x,f0]=fmincon_global(P,-10,10,5,10), toc % 求全局最优解 如果调用100 次这个求解程序，极有可能得出原问题的全局最优解。就本例而言，测 试的100 次全得出了全局最优解，总的测试时间大约38.7 s，平均每次0.387 s。 >> tic, X=[]; % 启动秒表，并设置空矩阵记录每次求解结果 for i=1:100 % 运行100 次求解函数，并评价找到全局最优解的成功率 [x,f0]=fmincon_global(P,-10,10,5,10); X=[X; x']; % 记录结果 end, toc % 显示100 次求解所需的总时间 3.4.4 最优曲线拟合方法 假设有一组数据xi, yi, i = 1, 2, · · · ,N，且已知这组数据满足某一函数原型 ?y(x) = f(a, x)，其中a待定系数向量，则最小二乘曲线拟合的目标就是求出这一组 待定系数的值，使得目标函数（即拟合的总误差） J = min a N Xi=1 yi ? ?y(xi)2 = min a N Xi=1 yi ? f(a, xi)2 （3-4-3） 为最小。在MATLAB 的最优化工具箱中提供了lsqcurvefit() 函数，可以解决最 小二乘曲线拟合的问题，该函数的调用格式为： [a,Jm]=lsqcurvefit(Fun,a0,x,y,xm,xM,options) ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. 其中，a0 为最优化的初值，x、y 为原始输入输出数据向量，Fun 为原型函数的 MATLAB 表示，可以用匿名函数描述，也可以用M-函数表示，该函数还允许指 定待定向量的最小值xm 和最大值xM，也可以设置搜索控制参数options。调用该 函数则将返回待定系数向量a，以及在此待定系数下的目标函数的值Jm。 例3-21 假设在实验中测出一组数据，且已知其可能满足的函数为 y(x) = a1ea2x sin(a3x + a4) cos(a5x) 其中ai 为待定系数，试用最小二乘方法拟合该函数。 解则可以通过最小二乘拟合的方法拟合出函数的待定系数。假设通过数据生成的方法 产生这组“实验数据”，下面将演示曲线的最小二乘拟合方法。首先由下面语句生成“实 验数据”，这些生成的坐标点可以用二维图形绘制出来，如图3-5（a）所示。 >> x=[0:0.01:0.1, 0.2:0.1:1,1.5:0.5:10]; % 生成不等间距的横坐标 y=0.56*exp(-0.2*x).*sin(0.8*x+0.4).*cos(-0.65*x); % 实验数据 plot(x,y,'o',x,y) % 绘制实验点坐标图形 0 2 4 6 8 10 --0.1 0 0.1 0.2 0.3 0.4 （a）给定的数据曲线 0 2 4 6 8 10 --0.1 0 0.1 0.2 0.3 0.4 （b）多项式拟合效果 图3-5 给定数据点的曲线拟合 由原型函数可以通过最小二乘进行最优拟合，这样可以通过MATLAB 语言编写匿 名函数，然后由最小二乘法求取待定系数。 >> F=@(a,x)a(1)*exp(-a(2)*x).*sin(a(3)*x+a(4)).*cos(-a(5)*x); f=optimset; f.RelX=1e-10; f.TolFun=1e-15; % 指定较高的拟合精度 a=lsqcurvefit(F,[1;1;1;1;1],x,y,[0,0,0,0,0],[],f) % 参数拟合 a0=[0.56;0.2;0.8;0.4;0.65]; norm(a-a0) % 和真值比较的误差 x0=0:0.01:10; y0=F(a0,x0); % 设置拟合点 y1=F(a,x0); plot(x0,y0,x0,y1,'--',x,y,'o') % 绘制拟合曲线 得出的拟合参数为a = [0.56, 0.2, 0.8, 0.4, 0.65]T，与给定的真值完全一致，拟合误差为 4.4177 × 10?7，拟合结果和图3-5（a）给出的完全一致。 MATLAB 提供的多项式拟合函数a=polyfit(x,y,n) 也可以用于曲线拟合，其 中x和y 为数据向量，n为拟合系数的阶次，通过该函数的调用将得出拟合多项式系数 向量a，该向量是多项式系数按降幂排列构成的向量，前面已经介绍过，用polyval() 函 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解101 数可以对多项式求值。下面的语句将比较6 次和8 次多项式拟合的效果，如图3-5（b）所 示，可见多项式拟合效果难以保证，故曲线拟合时如果已知原型时不必采用多项式拟合， 若不知原型时应该选择插值。 >> p=polyfit(x,y,6), y2=polyval(p,x0); % 6 次多项式拟合结果 p=polyfit(x,y,8); y3=polyval(p,x0); % 8 次多项式拟合 plot(x0,y0,x,y,'o',x0,y2,x0,y3) % 拟合结果比较 其中，拟合的6 次多项式拟合为 P6(x) = ?0.0002x6+0.0054x5?0.0632x4+0.3430x3 ? 0.8346x2+0.6621x+0.2017 3.5 Laplace变换与z 变换问题的MATLAB求解 在早期连续控制系统的研究中，常微分方程是最主要的建模工具，然而，由于 微分方程相对于代数方程要复杂得多，所以应该利用一种积分变换Laplace 变 换，将其映射成代数方程，从而引入了传递函数模型，该模型奠定了经典控制理论 的基础，线性系统时域响应解析解方法利用了Laplace 反变换的功能，而求解反变 换需要复变函数中的留数方法。同样，离散系统也可以利用z 变换构建离散传递函 数模型，求解传递函数又需要z 反变换。 3.5.1 Laplace变换 一个时域函数f(t) 的Laplace 变换可以定义为 L[f(t)] = Z∞ 0 f(t)e?stdt = F(s) （3-5-1） 其中L[f(t)] 为Laplace 变换的简单记号。 如果已知函数的Laplace 变换式F(s)，则可以通过下面的反变换公式直接求 出其Laplace 反变换 f(t) = L?1[F(s)] = 1 2πj Zσ+j∞ σ?j∞ F(s)estds （3-5-2） 其中σ 大于F(s) 奇异点的实部。 MATLAB 的符号运算工具箱提供的laplace() 函数及ilaplace() 函数直接 求取给定函数的正、反Laplace 变换。求解积分变换问题的步骤为： （1）声明符号变量。用syms 命令声明符号变量。 （2）描述原函数。直接用MATLAB 格式描述出原函数。 （3）求取积分变换。Laplace 变换和反变换可以由laplace() 和ilaplace() 直 接变换。其实，fourier()、ifourier() 函数可以求解Fourier 变换与反变换。 例3-22 试求时域函数f(t) = 1 ? (1 + at)e?at 的Laplace 变换。 解下面通过计算机工具直接求取这个函数的Laplace 变换。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. >> syms a t % 声明所需的变量为符号变量 f=1-(1+a*t)*exp(-a*t); % 表示时域函数公式 F=laplace(f), % 求取函数的Laplace 变换，注意得出的结果不一定最简 可以得出如下的结果 F = 1 s ? 1 s + a ? a (s + a)2 利用ilaplace() 对上述结果进行Laplace 反变换，则可以还原成原来的时域函数，可 以采用下面命令来完成这样的反变换。由ilaplace(F) 可以得出反变换的结果为1 ? e?at ? ate?at。 例3-23 已知某Laplace 表达式如下，试求Laplace 反变换。 G(s) = s + 3 s(s4 + 2s3 + 11s2 + 18s + 18) 解用下面的MATLAB 语句就可以求出信号Laplace 反变换的解析解。 >> syms s, G=(s+3)/s/(s^4+2*s^3+11*s^2+18*s+18); ilaplace(G) 可以得出解析解为 1 255 cos 3t ? 13 255 sin 3t ? 29 170 e?t cos t ? 3 170 e?t sin t + 1 6 例3-24 求出Laplace 变换式L?1  3a2 s3 + a3 , a > 0。 解如果想证明该式子，可以首先对给出的式子进行Laplace 反变换。 >> syms s t; syms a positive; F=3*a^2/(s^3+a^3); f=simplify(ilaplace(F)) 可以求出L?1  3a2 s3 + a3 = e?at + eat/2 ?cos √3 2 at + √3 sin √3 2 at。 例3-25 试求Laplace 函数F(s) = e?√s √s(√s ? 1) 的反变换。 解反变换可以通过下面函数直接求出 >> syms s; z=sqrt(s); f=ilaplace(exp(-z)/z/(z-1)) 在新版本下这个问题是不能求解的，如果使用早期版本，例如，MATLAB 2008a 及 更早的版本，则得出的结果为f(t) = et?1erfc??? (2t ? 1)/(2√t)，原函数是分数阶导 数的一个例子。 3.5.2 数值Laplace变换 前面给出了Laplace 变换的求解函数laplace()，该函数可以推导出很多时域 函数的Laplace 变换的解析解，同时也有很多函数Laplace 变换的解析解不存在或 不适合用解析解方法求解，所以应该考虑数值方法求解Laplace 变换问题。 Juraj Valsa 开发了基于数值方法的Laplace 反变换的函数INVLAP() [10,11]，该 函数的调用格式为[t,y]=INVLAP(f,t0,tn,N, 其他参数)，其中，原函数由含有字 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解103 符s 的字符串表示，(t0, tn) 为感兴趣的区间且t0 ?= 0，N 为用户选择的计算点数，用 户可以选择不同的N 值来检验运算的结果。“其他参数”的选取可以参考原函数的 联机帮助，不过这里建议：除非特别需要没有必要去人为修改这些默认参数。 文献[1]对INVLAP() 函数进行了扩展，给出了基于数值Laplace 变换的反馈控 制系统输出响应求解函数INVLAP_new()，其调用格式如下： [t,y]=INVLAP_new(G,t0,tn,N), %G的Laplace 反变换 [t,y]=INVLAP_new(G,t0,tn,N,H), %G，H 闭环系统的脉冲响应 [t,y]=INVLAP_new(G,t0,tn,N,H,u), % u用于描述输入信号 [t,y]=INVLAP_new(G,t0,tn,N,H,tx,ux), % tx，ux 为时间、输入采样点 该函数支持多种调用格式，其中，G为Laplace 变换表达式的字符串，如果同时 提供了H，则G为前向通路的传递函数模型字符串，H 为负反馈回路传递函数的字 符串；如果需要描述输入信号，则u可以为输入信号的Laplace 变换字符串，或输入 时域信号的匿名函数句柄；输入信号还可以由采样点(tx,ux) 表示；如果只考虑G 模型的响应，可以将H 设置为0。 例3-26 假设复杂开环无理模型如下[12]，试绘制单位负反馈统的阶跃响应曲线。 G(s) = sinh(0.1√s) 0.1√s 2 1 √s sinh(√s) 解开环无理传递函数可以由字符串表示，这样由下面语句直接绘制出系统的闭环阶跃 响应曲线，如图3-6 所示。 >> G='(sinh(0.1*sqrt(s))/0.1/sqrt(s))^2/sqrt(s)/sinh(sqrt(s))'; [t,y]=INVLAP_new(G,0,10,1000,1,'1/s'); plot(t,y) % 闭环系统阶跃响应 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 1.2 图3-6 闭环系统的阶跃响应曲线 3.5.3 z 变换 离散序列信号f(n) 的z 变换可以定义为 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. Z [f(n)] = ∞Xi=0 f(n)z?n = F(z) （3-5-3） 给定z 变换式子F(z)，则其z 反变换的数学表示为 f(n) = Z ?1[F(z)] = 1 2πj IF(z)zn?1dz （3-5-4） 可以调用MATLAB 符号运算工具箱中的ztrans() 函数和iztrans() 函数对 给定的函数进行正、反z 变换。 例3-27 一般介绍z 变换的书中不介绍q/(z?1 ? p)m（p ?= 0）函数的z 反变换，而该函 数是求取有重根离散系统解析解的基础。试推导该函数的z 反变换一般表达式。 解如果m选作符号变量，则ztrans() 函数并不能得出有益的结果。这里尝试对不同 的m值进行反变换，并总结出一般规律。根据要求，可以用符号运算工具箱求出m = 1, 2, · · · , 5 的z 反变换。 >> syms p q z n; assume(p~=0); assume(n~=0) for i=1:5 % 尝试不同阶次，以便总结规律 F=iztrans(q/(1/z-p)^i); % 求反变换 F=subs(F,{nchoosek(n-1,2),nchoosek(n-1,3),nchoosek(n-1,4)},... {(n-1)*(n-2)/2,(n-1)*(n-2)*(n-3)/6,... (n-1)*(n-2)*(n-3)*(n-4)/24}); F=prod(factor(F)) % 对得出的结果直接作因式分解 end 上述语句可以直接得出如下结果。 F1 = ? q p1+n , F2 = q p2+n (1 + n) , F3 = ? q 2p3+n (1 + n)(2 + n) F4 = q 6p4+n (3 + n) (2 + n) (1 + n) , F5 = ? q 24p5+n (4 + n) (3 + n) (2 + n) (1 + n) 总结上述结果的规律，可以写出一般的z 反变换结果 Z ?1  q (z?1 ? p)m= (?1)mq (m ? 1)! pn+m(n + 1)(n + 2) · · · (n + m ? 1) 注意，求解这个问题在新版本下使用的语句比较烦琐，因为，默认得出的结果含有 nchoosek() 函数，即组合函数，所以，需要手工替换成相应的多项式形式。此外，直接采 用simplify() 函数对结果化简，格式可能不统一，所以这里采用因式分解的方法化简 结果。早期版本中，对本问题只需给出命令simplify() 即可。 例3-28 已知某信号的z 变换表达式如下，试求其z 反变换。 H(z) = z(5z ? 2) (z ? 1)(z ? 1/2)3(z ? 1/3) 解可以通过下面的MATLAB 的符号运算工具箱直接求取该信号的z 反变换。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解105 >> syms z; H=z*(5*z-2)/((z-1)*(z-1/2)^3*(z-1/3)); iztrans(H) 得出的结果为 Z ?1[H(z)] = 36 + (72 ? 60n ? 12n2)1 2n ? 1081 3n 3.6 习题 （1）对下面给出的各个矩阵求取各种参数，如矩阵的行列式、迹、秩、特征多项式、范 数等，试分别求出它们的解析解。 A =264 7.5 3.5 0 0 8 33 4.1 0 0 9 103 ?1.5 0 0 3.7 19.3 375 , B =264 5 7 6 5 7 10 8 7 6 8 10 9 5 7 9 10 375 （2）求出下面给出的矩阵的秩和Moore–Penrose 广义逆矩阵，并验证它们是否满足 Moore–Penrose 逆矩阵的条件。 A =26664 2 2 3 1 2 2 3 1 4 4 6 2 1 1 1 1 ?1 ?1 ?1 3 37775 , B =24 4 1 2 0 1 1 5 15 3 1 3 535 （3）试求解下面的线性代数方程组，并检验解的正确性。 X264 7 6 9 7 7 1 3 2 2 1 5 5 6 4 2 6 375 = 2 1 0 1 0 3 1 2  （4）给定下面特殊矩阵A，试利用符号运算工具箱求出其逆矩阵、特征值，并求出状态 转移矩阵eAt 的解析解。 A =26664 ?9 11 ?21 63 ?252 70 ?69 141 ?421 1684 ?575 575 ?1149 3451 ?13801 3891 ?3891 7782 ?23345 93365 1024 ?1024 2048 ?6144 24573 37775 （5）试求下面齐次方程的基础解系。 8> >><>>>: 6x1 + x2 + 4x3 ? 7x4 ? 3x5 = 0 ?2x1 ? 7x2 ? 8x3 + 6x4 = 0 ?4x1 + 5x2 + x3 ? 6x4 + 8x5 = 0 ?34x1 + 36x2 + 9x3 ? 21x4 + 49x5 = 0 ?26x1 ? 12x2 ? 27x3 + 27x4 + 17x5 = 0 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. （6）求解下面的Lyapunov 方程，并检验所得解的精度。 24 1 2 3 4 5 6 7 8 0 35 X +X 24 2 3 6 3 5 2 3 2 2 35 = 24 1 3 2 3 4 1 5 2 1 35 （7）试用数值方法和解析方法求取下面的Sylvester 方程，并验证得出的结果。 26664 3 ?6 ?4 0 5 1 4 2 ?2 4 ?6 3 ?6 7 3 ?13 10 0 ?11 0 0 4 0 3 4 37775 X +X 24 3 ?2 1 ?2 ?9 2 ?2 ?1 9 35 = 26664 ?2 1 ?1 4 1 2 5 ?6 1 6 ?4 ?4 ?6 6 ?3 37775 （8）试求出下面的代数方程的全部的根 8> <> : x2y2 ? zxy ? 4x2yz2 = xz2 xy3 ? 2yz2 = 3x3z2 + 4xzy2 y2x ? 7xy2 + 3xz2 = x4zy （9）采用适当的方法求解下面的非线性方程[13] à 8> <> : xyz = 1 x2 + 2y2 + 4z2 = 7 2x2 + y3 + 6z = 7 á 8<: x2 + 2 sin(yπ/2) + z2 = 0 ?2xy + z = 3 x2z ? y = 7 （10）试找出下面非线性矩阵方程所有可能的解 24 9 1 0 8 7 3 3 0 6 35X2 +X 24 5 6 6 9 2 2 6 9 5 35 X +X 24 7 9 4 6 7 1 4 6 4 35 + I3×3 = 0 （11）若A、B、C 和D矩阵由例3-13 给出，试求解矩阵方程 X3 +X4D?X2BX + CX?I = 0 假设已经求出了方程的77 个实根，总共3351 个复数根，在data3ex1.mat 文件给 出，试接着求解该方程，看看能不能找到新的解。 （12）试求出伪多项式方程x √7 + 2x √3 + 3x √2?1 + 4 = 0 所有的根，并检验结果。 （13）某Riccati 方程的数学表达式为PA + A T P ? PBR?1B T P + Q = 0，且已知 A = 264 ?27 6 ?3 9 2 ?6 ?2 ?6 ?5 0 ?5 ?2 10 3 4 ?11 375 , B = 264 0 3 16 4 ?7 4 9 6 375 Q = 264 6 5 3 4 5 6 3 4 3 3 6 2 4 4 2 6 375 , R = 4 1 1 5  试求解该方程，得出P 矩阵，并检验得出解的精度。试求出并验证方程全部的根。 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解107 （14）试求解非线性矩阵方程AX sin(B2X + C)X +Xe?B + C = 0，其中 A =24 4 4 4 0 3 9 4 7 935 ,B =24 1 1 7 3 5 5 8 1 035 ,C =24 6 4 7 1 9 4 1 3 035 （15）Lorenz 方程是研究混沌问题的著名的非线性微分方程，其数学形式为 8> <> : x˙ 1(t) = ?βx1(t) + x2(t)x3(t) x˙ 2(t) = ?σx2(t) + σx3(t) x˙ 3(t) = ?x1(t)x2(t) + γx2(t) ? x3(t) 其中，β=8/3, σ=10, γ=28，且其初值为x1(0)=x2(0)=0, x3(0)=10?3。试求出 其数值解，绘制三维空间的相轨迹，并绘制出Lorenz 方程解在两两平面上的投影。 （16）请给出求解下面微分方程的MATLAB 命令 ... y + ty¨y + t2 y˙y2 = e?ty, y(0) = 2, y˙(0) = y¨(0) = 0 并绘制出y(t) 曲线，试问该方程存在解析解吗？ （17）Lotka–Volterra 扑食模型方程如下，且初值为x(0) = 2, y(0) = 3，试求解该微分 方程，并绘制相应的曲线。 ( x˙ (t) = 4x(t) ? 2x(t)y(t) y˙(t) = x(t)y(t) ? 3y(t) （18）试求解下面的零初值微分方程 à ( x˙ (t) = px2(t) ? y(t) + 3 ? 3 y˙(t) = arctan(x2(t) + 2x(t)y(t)) á ( x˙ (t) = ln(2 ? y(t) + 2y2(t)) y˙(t) = 4 ?px(t) + 4x2(t) （19）试求解边值问题¨y(x) = λ2(y2(x)+cos2 πx)+2π2 cos 2πx，其中，y(0) = y(1) = 0， y˙(0) = 1。提示：这个问题不是真正的边值问题，可以考虑由y(1) = 0 求出未知参 数λ，再求解微分方程。试学习与使用bvp5c() 函数直接求解这个边值问题。 （20）试求出下面微分方程组的解析解，并和数值解比较。 à ( x¨(t) = ?2x(t) ? 3x˙ (t) + e?5t, x(0) = 1, x˙ (0) = 2 y¨(t) = 2x(t) ? 3y(t) ? 4x˙ (t) ? 4y˙(t) ? sin t, y(0) = 3, y˙(0) = 4 á ( ¨x(t) + ¨y(t) + x(t) + y(t) = 0, x(0) = 2, y(0) = 1 2x¨(t) ? y¨(t) ? x(t) + y(t) = sin t, x˙ (0) = y˙(0) = ?1 （21）一级倒立摆模型的数学描述为 8>><> >: ¨x = u + ml sin θθ˙2 ? mg cos θ sin θ M + m ? mcos2 θ ¨θ = u cos θ ? (M + m)g sin θ + ml sin θ cos θθ˙ ml cos2 θ ? (M + m)l 已知m = M = 0.5 kg，l = 0.3 m，g = 9.81 m/s2。试求解该系统在单位阶跃信号u 作用下的零初始状态时间响应（注意：该系统为自然不稳定系统，如果需要使之稳 定则应使用特殊的控制方法）。 ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 控制系统计算机辅助设计——MATLAB 语言与应用（第4 版. .）. （22）假设方程为零初值问题，试求解下面的常微分方程 (cos ¨x(t) ... y (t) ? cos x¨(t) ? y¨(t) ? x(t)y˙(t) + e?x(t)y(t) = 2 sin ¨x(t) cos ... y (t) ? x(t)y˙(t) + x¨(t)y(t) ? y2(t)y˙(t) = 5 （23）求解下面的最优化问题。 à min x s.t.( 4x2 1+x2 2?4 x1,x2?0 x21 ? 2x1 + x2 á max x s.t.x1+x2+5=0h?(x1 ? 1)2 ? (x2 ? 1)2i （24）考虑下面二元最优化问题的求解，还可以用图解方法验证你得出的解。 max x s.t.( 9?x2 1+x2 2 x1+x2?1 (?x21 ? x2) （25）试求解下面的非线性规划问题。 min x s.t. 8> >>>><>>>>>: 0.003079x3 1x3 2x5?cos3 x6?0 0.1017x3 3x3 4?x2 5 cos3 x6?0 0.09939(1+x5)x3 1x2 2?cos2 x6?0 0.1076(31.5+x5)x3 3x2 4?x2 5 cos2 x6?0 x3x4(x5+31.5)?x5[2(x1+5) cos x6+x1x2x5]?0 0.2?x1?0.5,14?x2?22,0.35?x3?0.6 16?x4?22,5.8?x5?6.5,0.14?x6?0.2618 1 2 cos x6 x1x2(1 + x5) + x3x4 1 + 31.5 x5  （26）试求解下面的最优化问题[14]。 min q,w,k s.t. 8> >><>>>: q3+9.625q1w+16q2w+16w2+12?4q1?q2?78w=0 16q1w+44?19q1?8q2?q3?24w=0 2.25?0.25k?q1?2.25+0.25k 1.5?0.5k?q2?1.5+0.5k 1.5?1.5k?q3?1.5+1.5k k （27）试求解下面的最优化问题。 min x s.t. 8> >><>>>: 0.0193x3?x1?0 0.00954x3?x2?0 750×1728?πx2 3x4?4πx3 3/3?0 x4?240?0 0.0625?x1,x2?6.1875,10?x3,x4?200 0.6224x1x2x3x4 + 1.7781x2x23 + 3.1661x21 x4 + 19.84x21 x3 （28）试用一般的有约束最优化方法求解下面的线性规划问题，试学习linprog() 函数 的使用方法重新求解下面的问题。 à min x s.t.8>><>>: 4x1?x2+2x3?x4=?2 x1+x2?x3+2x4?14 2x1?3x2?x3?x4??2 x1,2,3??1,x4无约束 ?3x1 + 4x2 ? 2x3 + 5x4 á min x s.t.8>><>>: x1+x2+x3+x4=4 ?2x1+x2?x3?x6+x7=1 3x2+x3+x5+x7=9 x1,2,··· ,7?0 x6 + x7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .? 第3 章科学运算问题的MATLAB 求解109 （29）考虑一个简单的一元函数最优化问题求解，f(x) = x sin(10πx) + 2，x ∈ (?1, 2)， 试求出f(x) 取最大值时x 的值。已知，该函数曲线有很强振荡，所以采用常规 最优化方法时，若初值选择不当往往会得出局部最小值。要求在本题求解中，在 x ∈ (?1, 2) 区间内随机选择40 个初始点，按照图3-7 中给出的流程编程，用循环 的方式从每个初始点出发进行搜索，得出全局最优解。如有可能，将该方法和函数 扩展成求取多元函数y = f(x) 全局最优解的通用程序。 选择40 个初始点向量x 令ym = 100 取x0 为x 向量值循环 调用fminsearch() 函数由初 值x0 求最优解，返回x1, y1 循环结束，输出全局最 优解x 和最小值ym 比较y1 和ym 大小，若y1