北航最优化方法大作业参考
1 流量工程问题
1.1 问题重述
定义一个有向网络G=(N,E)，其中N是节点集，E是弧集。令A是网络G的点弧关联矩阵，即N×E阶矩阵，且第l列与弧里(I,j)对应，仅第i行元素为1，第j行元素为-1，其余元素为0。再令b m=(b m1,…,b mN)T，f m=(f m1,…,f mE)T，则可将等式约束表示成：
Af m=b m
本算例为一经典TE算例。算例网络有7个节点和13条弧，每条弧的容量是5个单位。此外有四个需求量均为4个单位的源一目的对，具体的源节点、目的节点信息如图所示。这里为了简单，省区了未用到的弧。此外，弧上的数字表示弧的编号。此时，c=((5,5…,5)1 )T，
×13
根据上述四个约束条件，分别求得四个情况下的最优决策变量x=((x12,x13,…,x75)1×13)。
图 1 网络拓扑和流量需求
1.2 7节点算例求解
1.2.1 算例1（b1=[4;-4;0;0;0;0;0]T）
转化为线性规划问题：
Minimize c T x1
Subject to Ax1=b1
x1>=0 利用Matlab编写对偶单纯形法程序，可求得:
最优解为x1*=[4 0 0 0 0 0 0 0 0 0 0 0 0]T
对应的最优值c T x1=20
1.2.2 算例2（b2=[4;0;-4;0;0;0;0]T）
Minimize c T x2
Subject to Ax2=b2
X2>=0 利用Matlab编写对偶单纯形法程序，可求得:
最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T
对应的最优值c T x2=20
1.2.3 算例3（b3=[0;-4;4;0;0;0;0]T）
Minimize c T x3
Subject to Ax3=b3
X3>=0 利用Matlab编写对偶单纯形法程序，可求得:
最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T
对应的最优值c T x3=40
1.2.4 算例4（b4=[4;0;0;0;0;0;-4]T）
Minimize c T x4
Subject to Ax4=b4
X4>=0
利用Matlab编写对偶单纯形法程序，可求得:
最优解为x4*=[4 0 0 4 0 0 0 0 0 4 0 0 0]T
对应的最优值c T x4=60
1.3 计算结果及结果说明
1.3.1 算例1（b1=[4;-4;0;0;0;0;0]T）
算例1中，由b1可知，节点2为需求节点，节点1为供给节点，由节点1将信息传输至节点2的最短路径为弧1。
图 2 算例1最优传输示意图
求得的最优解为x1*=[4 0 0 0 0 0 0 0 0 0 0 0 0]T，即只经过弧1运输4个单位流量，其余弧无流量。又因为，每条弧的费用均为5，所以最小费用为20。
经分析，计算结果合理可信。
1.3.2 算例2（b2=[4;0;-4;0;0;0;0]T）
算例2中，由b2可知，节点3为需求节点，节点1为供给节点，由节点1将信息传输至节点2的最短路径为弧2。
图 3 算例2最优传输示意图
求得的最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T，即只经过弧2运输4个单位流量，其余弧无流量。又因为，每条弧的费用均为5，所以最小费用为20。
经分析，计算结果合理可信。
1.3.3 算例3（b3=[0;-4;4;0;0;0;0]T）
算例3中，由b3可知，节点2为需求节点，节点3为供给节点，由节点3将信息传输至节点2的最短路径为弧5->弧1。
图 4 算例3最优传输示意图
求得的最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T，即经过弧5运输4个单位流量至节点1，再经弧1运输4个单位流量至节点2，其余弧无流量。又因为，每条弧的费用均为5，所以最小费用为40。
经分析，计算结果合理可信。
1.3.4 算例4（b4=[4;0;0;0;0;0;-4]T）
算例4中，由b4可知，节点7为需求节点，节点1为供给节点，由节点1将信息传输至节点7的最短路径为弧1->弧4->弧10。
图 5 算例3最优传输示意图
求得的最优解为x4*=[4 0 0 4 0 0 0 0 0 4 0 0 0]T，即经过弧1运输4个单位流量至节点2，再经弧4运输4个单位流量至节点5，最后经弧5运输4个单位流量至节点7，其余弧无流量。又因为，每条弧的费用均为5，所以最小费用为60。
经分析，计算结果合理可信。
2 重要算法编写与观察
2.1 习题5.6
（a）初值为（0,0）时
本算法令g的2范数在<10-4时，停止迭代，经过86次迭代收敛。
收敛因子（f(k+1)-f*）/（f(k)-f*）=0.7623
图 6 收敛因子截图
（b）初值为（-0.4,0）时
本算法令g的2范数在<10-4时，停止迭代，经过112次迭代收敛。
收敛因子（f(k+1)-f*）/（f(k)-f*）=0.81
图7 收敛因子截图
（c）初值为（10,0）时
本算法令g的2范数在<10-4时，停止迭代，经过5次迭代收敛。
收敛因子（f(k+1)-f*）/（f(k)-f*）=3.9022e-4
图8 收敛因子截图
（d）初值为（11,0）时
本算法令g的2范数在<10-4时，停止迭代，经过2次迭代收敛。
收敛因子（f(k+1)-f*）/（f(k)-f*）= 0
图9 收敛因子截图
图10 自变量（x1,x2）截图
总结：最速降线法的收敛因子随着初值的不同而变化，对于个别初值（如本习题初值取（11,0）时），算法可迅速收敛。因此，初值的选取对于最速降线法的收敛速度有较大影响。 2.2 习题5.7
（a ） 由()94ln(7)f x x x =--可得：
4
'()97
f x x =-
- 2
4
"()(7)
f x x =
- 故，牛顿迭代法的确切公式为：
2
4974(7)x s x -
-=-
- （b ）从以下五个初值开始迭代 （1）x(0)=7.40
（2）
x(0)=7.20
（3）x(0)=7.01
（4）x(0)=7.80
（5）x(0)=7.88
（c）本问题的最优值为7.4444444。由上述五个初值点的前五步迭代可以看出：当初值点在区间（7.4444444，7.8888）内时，第二次迭代点将落在（7，7.4444444）之间，随后逐渐增加，直至逼近最优值。
当初值点在区间（7，7.4444444）内时，则迭代点逐渐增加，逼近最优值。
当取初值不在（7,7.8888）内时，牛顿法不收敛。
2.3 习题5.8
（a）没有线搜索的牛顿法
μ=0.1时，
μ=1时，
（b）具有线搜索的牛顿法
μ=0.1时，
μ=1时，
（未完成）
2.4 习题5.9
（a）初值选（1.2,1.2）时，
◆最速降线法：
本算法令g的2范数在<10-2时，停止迭代，经过3262次迭代得到以下结果。
图11 最速降线法初值为（1.2,1.2）的等值线图及迭代轨迹
◆牛顿法：
本算法令s的4范数在<10-6时，停止迭代，经过4次迭代得到以下结果。
图12 牛顿法初值为（1.2,1.2）的等值线图及迭代轨迹
（b）初值选（-1.2,1）时，
最速降线法：
本算法令g的4范数在<10-2时，停止迭代，经过6835次迭代得到以下结果。
图13 最速降线法初值为（-1.2,1）的等值线图及迭代轨迹
牛顿法：
本算法令s的4范数在<10-6时，停止迭代，经过6次迭代得到以下结果。
图14 牛顿法初值为（-1.2,1）的等值线图及迭代轨迹
2.5 习题5.19
N=5
迭代6次后，满足收敛条件。
N=8
迭代19次后，满足收敛条件。
N=14
迭代49次后，满足收敛条件。（表略）
N=40
迭代74次后，满足收敛条件。（表略）
2.6 习题5.27
调用MATLAB自带的lsqnonlin.m函数，计算可得对应的x(1)、x(2)和标准差如下表所示。
由上可知，标准差值较为恒定，随初值变化不十分显著；x1和x2值随初值选取的不同而不同。
2.7 习题6.4
（未完成）
3 附录
3.1 对偶单纯形法函数MATLAB程序
function [sol,val,kk]=duioudanchun(A,N)
B=A;
[mA,nA]=size(A);
kk=0;
flag=1;
while flag
kk=kk+1;
if A(:,nA)>=0
flag=0;
sol=zeros(1,nA);
for i=1:mA-1
sol(N(i))=A(i,nA);
end
val=sol*(B(mA,:))';
else
for i=1:mA-1
if A(i,nA)<0&A(i,1:nA-1)>=0
disp('have infinite solution!');
flag=0;
break;
end
end
if flag
temp=0;
for i=1:mA-1
if A(i,nA) temp=A(i,nA); outb=i; end end sita=zeros(1,nA-1); for i=1:nA-1 if A(outb,i)<0 sita(i)=A(mA,i)/A(outb,i); end end temp=-inf; for i=1:nA-1 if sita(i)<0&sita(i)>temp temp=sita(i); inb=i; end end for i=1:mA-1 if i==outb N(i)=inb; end end A(outb,:)=A(outb,:)/A(outb,inb); for i=1:mA if i~=outb A(i,:)=A(i,:)-A(outb,:)*A(i,inb); A(mA,nA)=0; end end end end end 3.2 最速降线法求Rosenbrock函数最小值matlab程序如下： function rb = rbfun(x,y) rb=100*(y-x^2)^2+(1-x)^2 end clear clc syms x y g G g=gradient(rb(x,y),[x y]) %定义梯度向量 G=hessian(rb(x,y),[x y]) %定义海森阵 X(1,:)=[-1.4 1]; %定义初始点 x=X(1,1);y=X(1,4); A(1,:)=subs(g) %给梯度赋初值 i=1 while(norm(A(i,:),4)>10^(-4)) %收敛条件 f(i)=rb(x,y) %记录函数值 P(i,:)=-A(i,:) %得到迭代方向 fz(i)=-A(i,:)*P(i,:)' %-gT*p %精确搜索法步长的分子fm(i)=P(i,:)*subs(G)*P(i,:)' %精确搜索法步长的分母a(i)=fz(i)/fm(i) %精确搜索法步长 X(i+1,:) = X(i,:)+a(i)*P(i,:) %产生新的点 x=X(i+1,1);y=X(i+1,4) A(i+1,:)=subs(g) %产生新的梯度 i=i+1 end 3.3 牛顿法求Rosenbrock函数最小值matlab程序如下： function rb = rbfun(x,y) rb=100*(y-x^2)^2+(1-x)^2 end clear clc syms x y g G g=gradient(rb(x,y),[x y]) %定义梯度向量 G=hessian(rb(x,y),[x y]) %定义海森阵 X(1,:)=[-1.4 1]; %定义初值 x=X(1,1);y=X(1,4); A(1,:)=subs(g) %给梯度赋初值 H=subs(inv(G)) %得到海森阵初值 S(1,:)=-A(1,:)*H %得到s初值 i=1 while (norm(S(i,:),4)>10^(-6)) %收敛条件 f(i)=rb(x,y) %定义函数值 X(i+1,:) = X(i,:)+S(i,:) %得到下一迭代点x=X(i+1,1);y=X(i+1,4) %给x,y分别赋值 A(i+1,:)=subs(g) %得到新的梯度值H=subs(inv(G)) %得到新的海森阵S(i+1,:)=-A(i+1,:)*H %得到新的增量s i=i+1 end 3.4 共轭梯度法求解习题5.19程序如下： clear clc K=40 G=zeros(K,K) for m = 1: K for n = 1: K G(m,n)=1/(m+n-1) end end X(1,:)=zeros(1,K) b=ones(1,K) A(1,:)=X(1,:)*G-b P(1,:)=-A(1,:) i=1 while (norm(A(i,:),4)>10^(-6)) %收敛条件 d=P(i,:)*G fz(i)=A(i,:)*A(i,:)' %精确搜索法步长的分子 fm(i)=P(i,:)*d' %精确搜索法步长的分母 a(i)=fz(i)/fm(i) %精确搜索法步长X(i+1,:) = X(i,:)+a(i)*P(i,:) %产生新的点 A(i+1,:)=A(i,:)+a(i)*d beta(i+1)=(A(i+1,:)*A(i+1,:)')/(A(i,:)*A(i,:)') P(i+1,:)=-A(i+1,:)+beta(i+1)*P(i,:) i=i+1 end ①图与网 破圈法：任取一个圈，去掉一条权最大的边，直到最小树。 避圈法：选最小权的边，避圈前进，直到最小树。 最短路算法： Dijkstra法：从V s给定P标号T标号λ标号（T标号变为P标号λ标号记位置） 反向追踪：列表，d1(V1,V j)→d k(V1,V j)=min(ωij+d k(V1,V i))据最小权反向追踪 网络优化： 最小截集最大流：找到最小截集（弧的集合） 标号法：开始，为的标号， 最小费用最大流： 邮递员问题：通过消灭奇点，找欧拉回路 网络计划图： 最早开始最晚开始机动时间 最早结束最晚结束自由时差 工期优化：人力，费用，工期优化。 费用率=（最短时间费用-正常时间费用）/（正常时间-最短时间）②排队论（保证服务质量，又减少费用） 顾客源→（排队规则）队列→（服务规则）服务机构→离去 服务规则：FCFS，LCFS，随机服务，PR M(顾客到达)|A(服务时间)|1(服务台数)|∞（容量）|∞（顾客源） N(t)队长N q (t)排队长T(t)顾客逗留时间T q (t)顾客等待时间 L 平均队长L q 平均等待队长W 平均逗留时间W q 平均等待时间 R 为系统利用率 泊松流(M)：无后效性；平稳性；单个性； P 1(t,t+Δt)=λΔt+o(Δt)； o(Δt)=∑∞ 2P n (t,t+Δt);E ξ=D ξ=λt （t 时刻n 个顾客的概率） 负指数分布(M)：无记忆性(P(T>t+s/t>s)=P(T>t));[0,t)至少到达一 个顾客1-P 0(t ）=1-e -t λ,t>0 !)()(K t e t V K t k λλ-= ,2,1,0=K ?? ?<≥-=-0,00,1)(t t e t F t i λξ),2,1( =i 爱尔朗分布(E K )：（相当于泊松流到达后被k 个服务台均分顾客形成） (其中，t>0,E(T)=1/μ,Var(T)=1/μ2k ） )! 1()()(1 >-= --t e k t t f t k μμμ K=1为M ，k=∞定长分布D,k ≥30正态分布近似 G 表示一般相互独立的随机分布 Little 公式：（四者知一即可） μ1 + =q W W W L λ= q q W L λ= ρ+=q L L ∑∞ ==0 n n nP L ∑∑∞=∞ =+=-=s n n m s n q nP P s n L 0 )( 服务率：ρ=λ/μ(λ为到达μ为服务) 排队系统分析： 北航最优化方法大作业参考 1 流量工程问题 1.1 问题重述 定义一个有向网络G=(N,E)，其中N是节点集，E是弧集。令A是网络G的点弧关联矩阵，即N×E阶矩阵，且第l列与弧里(I,j)对应，仅第i行元素为1，第j行元素为-1，其余元素为0。再令b m=(b m1,…,b mN)T，f m=(f m1,…,f mE)T，则可将等式约束表示成： Af m=b m 本算例为一经典TE算例。算例网络有7个节点和13条弧，每条弧的容量是5个单位。此外有四个需求量均为4个单位的源一目的对，具体的源节点、目的节点信息如图所示。这里为了简单，省区了未用到的弧。此外，弧上的数字表示弧的编号。此时，c=((5,5…,5)1 )T， ×13 根据上述四个约束条件，分别求得四个情况下的最优决策变量x=((x12,x13,…,x75)1× )。 13 图 1 网络拓扑和流量需求 1.2 7节点算例求解 1.2.1 算例1（b1=[4;-4;0;0;0;0;0]T） 转化为线性规划问题： Minimize c T x1 Subject to Ax1=b1 x1>=0 利用Matlab编写对偶单纯形法程序，可求得: 最优解为x1*=[4 0 0 0 0 0 0 0 0 0 0 0 0]T 对应的最优值c T x1=20 1.2.2 算例2（b2=[4;0;-4;0;0;0;0]T） Minimize c T x2 Subject to Ax2=b2 X2>=0 利用Matlab编写对偶单纯形法程序，可求得: 最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T 对应的最优值c T x2=20 1.2.3 算例3（b3=[0;-4;4;0;0;0;0]T） Minimize c T x3 Subject to Ax3=b3 X3>=0 利用Matlab编写对偶单纯形法程序，可求得: 最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T 对应的最优值c T x3=40 1.用薄钢板制造一体积5m 3，长度不小于4m ，无上盖的货箱，要求钢板耗量最小。确定货箱的长x 1、宽x 2和高x 3。试列出问题的数学模型。 解：min 32312122x x x x x x z ++= s.t 5321=x x x 41≥x 0,,321≥x x x 2.将下面的线性规划问题表示为标准型并用单纯形法求解 max f=x 1+2x 2+x 3 s ．t ．2x 1+x 2-x 3≤2 -2x 1+x 2-5x 3≥-6 4x 1+x 2+x 3≤6 x i ≥0 i=1，2，3 解：先化标准形： Min 321x x x z -+= 224321=+-+x x x x 6525321=++-x x x x 646321=+++x x x x 列成表格： 1 2 1 610011460105122001112----- 可见此表已具备1°，2°，3°三个特点，可采用单纯形法。首先从底行中选元素-1，由2/2,6/2,6/4最小者决定选第一行第一列的元素2，标以记号，迭代一次得 1 2 1 2102310401162010021212 11-------- 再从底行中选元素-2/3，和第二列正元素1/2，迭代一次得 1 2 12 32 30 210231040116201002121211- ------ 再从底行中选元素-3，和第二列正元素2，迭代一次得 4 2 3 3 410120280114042001112--- 再迭代一次得 10 2 30 2 10 6 221023 1010213000421021013-- 选取最优解： 陆吾生教授是加拿大维多利亚大学电气与计算机工程系 (Dept. of Elect. and Comp. Eng. University of Victoria) 的正教授, 且为我校兼职教授，曾多次来我校数学系电子系讲学。陆吾生教授的研究方向是：最优化理论和小波理论及其在1维和2维的数字信号处理、数字图像处理、控制系统优化方面的应用。 现陆吾生教授计划在 2007 年 10－11 月来校开设一门为期一个月的短期课程“最优化理论及其应用”（每周两次，每次两节课），对象是数学系、计算机系、电子系的教师、高年级本科生及研究生，以他在2006年出版的最优化理论的专著作为教材。欢迎数学系、计算机系、电子系的研究生及高年级本科生选修该短期课程，修毕的研究生及本科生可给学分。 上课地点及时间：每周二及周四下午2:00开始，在闵行新校区第三教学楼326教室。（自10月11日至11月8日） 下面是此课程的内容介绍。 ----------------------------------- 最优化方法及应用 I. 函数的最优化及应用 1.1 无约束和有约束的函数优化问题 1.2 有约束优化问题的Karush-Kuhn-Tucker条件 1.3 凸集、凸函数和凸规划 1.4 Wolfe对偶 1.5 线性规划与二次规划 1.6 半正定规划 1.7 二次凸锥规划 1.8 多项式规划 1.9解最优化问题的计算机软件 II 泛函的最优化及应用 2.1 有界变差函数 2.2 泛函的变分与泛函的极值问题 2.3 Euler-Lagrange方程 2.4 二维图像的Osher模型 2.5 泛函最优化方法在图像处理中的应用 2.5.1 噪声的消减 2.5.2 De-Blurring 2.5.3 Segmentation ----------------------------------------------- 注：这是一门约二十学时左右的短期课程，旨在介绍函数及泛函的最优化理论和方法，及其在信息处理中的应用。只要学过一元及多元微积分和线性代数的学生就能修读并听懂本课程。课程中涉及到的算法实现和应用举例都使用数学软件MATLAB 华东师大数学系 发动机空燃比控制器 引言：我主要从事自动化相关研究。这里介绍我曾经接触过的发动机空燃比控制器设计中的优化问题。 发动机空燃比控制器设计中的最优化问题 AFR =a f m m && (1) 空燃比由方程（1）定义，在发动机运行过程中如果控制AFR 稳定在14.7可以获 得最好的动力性能和排放性能。如果假设进入气缸的空气流量a m &可以由相关单元检测得到，则可以通过控制进入气缸的燃油流量f m &来实现空燃比的精确控制。由于实际发动机的燃油喷嘴并不是直接对气缸喷燃油，而是通过进气歧管喷燃油，这么做会在进 气歧管壁上液化形成油膜，因此不仅是喷嘴喷出的未液化部分燃油会进入气缸，油膜 蒸发部分燃油也会进入气缸，如方程（2）。这样如何更好的喷射燃油成为了一个问题。 1110101122211ττττ?? ?? -?? ??????????=+????????-????????????-???? ? ??? ?? ????????? ?f f f v X x x u x x X x y =x && (2) 其中12、,==ff fv x m x m &&=f y m &,=fi u m &这里面，表示油膜蒸发量ff m &、fv m &表示为液化部分燃油、fi m &表示喷嘴喷射的燃油，在τf 、τv 、X 都已知的情况下，由现代控制理论知识，根据系统的增广状态空间模型方程（3） 0000001 1 011011114.70ττττ????-?? ??????????=-+-??????????????? ??????????????? ?? ??=?????? f f v v a X X u +q q m y q x x x &&& (3) 其中()0 14.7?t a q = y -m &。由极点配置方法，只要设计控制器方程（4），就可以 使得y 无差的跟踪阶跃输入，那么y 也能较好的跟踪AFR *a m /&。 12-- u =K q K x (4) 这里面的12、K K 确定，可由主导极点概念降维成两个参数12C ,C ，虽然都是最终稳态无差，但是目标是使得瞬态过程中y 和阶跃输入y r 的差异尽可能的小。所以原问 北航计算机控制系统大作业 ————————————————————————————————作者：————————————————————————————————日期: 计算机控制系统 大作业 姓名：陈启航 学号： 教师:周锐 日期:2016年6月1日 综合习题1 已知: 4 4)(+= s s D , 1） 试用 Z 变换、一阶向后差分、向前差分、零极点匹配、Tus ｔｉn 变换和预修正的Tus ｔin （设关键频率＝4）变换等方法将D （s)离散化，采样周期分别取为０.１s 和 0.4s; ２） 将 Ｄ（z ）的零极点标在Z 平面图上 3） 计算D （j ω)和各个Ｄ(e j ωT )的幅频和相频特性并绘图,ｗ由０~ 20r ad ，计算40 个点,应包括=4 点,每个T 绘一张图（Z 变换方法单画) 4） 计算 D(s)及T=0．1,T=０.4 时D （z ）的单位脉冲响应，运行时间为4 秒 5) 结合所得的结果讨论分析各种离散化方法的特点 6) 写出报告，附上结果。 解： （1） Z 变换法: a.离散化： T =0.1s 时， D (z )= 4z z ?0.6703； T =0.4s 时, D (z )= 4z z ?0.2019 ； ｂ.D (z )的零极点 c. D (jω)和D(e jωT )幅频相频特性曲线 连续系统: -1 -0.8-0.6-0.4-0.200.20.40.60.81 -1-0.8-0.6-0.4-0.200.20.40.60.81 零点 T=0.1s 时极点T=0.4s 时极点 T=0.1s时 T=0.4s时 惯性导航基础课程大作业报告（一）光纤陀螺误差建模与分析 班级：111514 姓名： 学号 2014年5月26日 一.系统误差原理图 二.系统误差的分析 （一）漂移引起的系统误差 1. εx ，εy ，εz 对东向速度误差δVx 的影响 clc;clear all; t=1:0.01:25; g=9.8; L=pi/180*39; Ws=2*pi/84.4*60; Wie=2*pi/24; R=g/(Ws)^2; e=0.1*180/pi; mcVx1=e*g*sin(L)/(Ws^2-Wie^2)*(sin(Wie*t)-Wie*sin(Ws*t)/Ws); mcVx2=e*((Ws^2-(Wie^2)*((cos(L))^2))/(Ws^2-Wie^2)*cos(Ws*t)-(Ws^2)*((sin(L))^2)*cos(Wi e*t)/(Ws^2-Wie^2)-(cos(L))^2); mcVx3=(sin(L))*(cos(L))*R*e*((Ws^2)*cos(Wie*t)/(Ws^2-Wie^2)-(Wie^2)*cos(Ws*t)/(Ws^2-Wi e^2)-1); plot(t,[mcVx1',mcVx2',mcVx3']); title('Ex,Ey,Ez 对Vx 的影响'); xlabel('时间t'); ylabel('Vx(t)'); 0,δλδL ,v v δδ legend('Ex-mcVx1','Ey-mcVx2','Ez-mcVx3'); grid; axis square; 分析：εx，εy，εz对东向速度误差δVx均有地球自转周期的影响，εx，εy还会有舒勒周期分量的影响，其中，εy对δVx的影响较大。 2.εx，εy，εz对东向速度误差δVy的影响 clc;clear all; t=1:0.01:25; g=9.8; L=pi/180*39; Ws=2*pi/84.4*60; Wie=2*pi/24; R=g/(Ws)^2; e=0.1*180/pi; mcVy1=e*g*(cos(Wie*t)-cos(Ws*t))/(Ws^2-Wie^2); mcVy2=g*sin(L)*e/(Ws^2-Wie^2)*(sin(Wie*t)-Wie/Ws*sin(Ws*t)); mcVy3=g*cos(L)*e/(Ws^2-Wie^2)*(sin(Wie*t)-Wie/Ws*sin(Ws*t)); plot(t,[mcVy1',mcVy2',mcVy3']); title('Ex,Ey,Ez对Vy的影响'); xlabel('时间t'); ylabel('Vy(t)'); legend('Ex-mcVy1','Ey-mcVy2','Ez-mcVy3'); grid; axis square; 航空科学与工程学院 《航空工程大型通用软件应用》大作业 机翼结构设计与分析 组号第3组 小组成员11051090 赵雅甜 11051093 廉佳 11051100 王守财 11051108 刘哲 11051135 张雄健 11051136 姜南 6月 目录 一 CATIA部分....................................... 错误!未定义书签。( 一) 作业要求..................................... 错误!未定义书签。( 二) 作业报告..................................... 错误!未定义书签。 1、三维模型图................................... 错误!未定义书签。 2、工程图....................................... 错误!未定义书签。 二 FLUENT部分...................................... 错误!未定义书签。( 一) 作业要求..................................... 错误!未定义书签。( 二) 作业报告..................................... 错误!未定义书签。 1、计算方法和流程............................... 错误!未定义书签。 2、网格分布图................................... 错误!未定义书签。 3、气动力系数................................... 错误!未定义书签。 4、翼型表面压力曲线............................. 错误!未定义书签。 5、翼型周围压力云图............................. 错误!未定义书签。 6、翼型周围x方向速度云图....................... 错误!未定义书签。 7、翼型周围y方向速度云图....................... 错误!未定义书签。 8、翼型周围x方向速度矢量图..................... 错误!未定义书签。 9、翼型周围y方向速度矢量图..................... 错误!未定义书签。 10、流线图...................................... 错误!未定义书签。 三 ANSYS部分....................................... 错误!未定义书签。( 一) 作业要求..................................... 错误!未定义书签。( 二) 作业报告..................................... 错误!未定义书签。 1、机翼按第一强度理论计算的应力云图............. 错误!未定义书签。 2、机翼按第二强度理论计算的应力云图............. 错误!未定义书签。 3、机翼按第三强度理论计算的应力云图............. 错误!未定义书签。 4、机翼按第四强度理论计算的应力云图............. 错误!未定义书签。 计算机控制系统 大作业 姓名：陈启航 学号： 教师：周锐 日期：2016年6月1日 综合习题1 已知: 4 4 )(+= s s D ， 1） 试用 Z 变换、一阶向后差分、向前差分、零极点匹配、Tustin 变换和 预修正的Tustin （设关键频率＝4）变换等方法将D(s)离散化，采样周期分别取为 和 ； 2） 将 D(z)的零极点标在Z 平面图上 3） 计算D （j ω）和各个D(e j ωT )的幅频和相频特性并绘图,w 由0~ 20ra d ，计算40 个点，应包括＝4 点，每个T 绘一张图（Z 变换方法单画） 4） 计算 D(s)及T=，T= 时D(z)的单位脉冲响应，运行时间为4 秒 5） 结合所得的结果讨论分析各种离散化方法的特点 6） 写出报告，附上结果。 解： （1） Z 变换法： a.离散化： T =0.1s 时， D (z )=4z z ?0.6703 ； T =0.4s 时， D (z )=4z z ?0.2019 ； b.D (z )的零极点 c. D (jω)和D(e jωT )幅频相频特性曲线 连续系统： T =0.1s 时 T =0.4s 时 d. D(s)和D(z)单位脉冲响应 D(s)单位脉冲响应： D(z)单位脉冲响应： T=0.1s时 T=0.4s时 （2）各种离散化方法： a.离散化后的D(z) 1、一阶向后差分： T=0.1s时 D(z)= 0.2857z z?0.7143 T=0.4s时 D(z)= 0.6154z z?0.3846 2、一阶向前差分：T=0.1s时 D(z)= 0.4 z?0.6 T=0.4s时 D(z)= 1.6 z+0.6 3、零极点匹配T=0.1s时 D(z)=0.1648(z+1) z?0.6703 T=0.4s时 D(z)=0.3991(z+1) z?0.2019 4、Tustin变换T=0.1s时 D(z)=0.1667(z+1) z?0.6667 T=0.4s时 D(z)= 0.4444(z+1) 5、预修正的Tustin变换（设关键频率=4） T=0.1s时 D(z)=0.1685(z+1) z?0.6629 T=0.4s时 D(z)=0.5073(z+1) z+0.0146 b.D(z)的零极点 1、一阶向后差分 习题二包括题目：P36页5（1）（4） 5（4） 习题三 包括题目：P61页1(1)(2); 3; 5; 6; 14；15(1) 1(1)(2)的解如下 3题的解如下 5,6题 14题解如下 14. 设22121212()(6)(233)f x x x x x x x =+++---, 求点在(4,6)T -处的牛顿方向。 解：已知 (1) (4,6)T x =-，由题意得 121212212121212(6)2(233)(3)()2(6)2(233)(3)x x x x x x x f x x x x x x x x +++-----?? ?= ?+++-----?? ∴ (1)1344()56g f x -?? =?= ??? 21212122211212122(3)22(3)(3)2(233)()22(3)(3)2(233)22(3)x x x x x x x f x x x x x x x x +--+--------? ??= ? +--------+--?? ∴ (1)2(1)1656()()564G x f x --?? =?= ?-?? (1)1 1/8007/400()7/4001/200G x --?? = ?--?? ∴ (1)(1)11141/100()574/100d G x g -?? =-= ?-?? 15（1）解如下 15. 用DFP 方法求下列问题的极小点 （1）22 121212min 353x x x x x x ++++ 解：取 (0) (1,1)T x =，0H I =时，DFP 法的第一步与最速下降法相同 2112352()156x x f x x x ++???= ?++??， (0)(1,1)T x =，(0) 10()12f x ???= ??? (1)0.07800.2936x -??= ?-??， (1) 1.3760() 1.1516f x ???= ?-?? 以下作第二次迭代 (1)(0) 1 1.07801.2936x x δ-??=-= ?-??， (1)(0) 18.6240()()13.1516f x f x γ-??=?-?= ?-?? 0110 111011101 T T T T H H H H H γγδδδγγγ=+- 目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。 这样即求出所有A所有实特征值对应的一个特征向量。 #include x zD 天津大学《最优化方法》复习题(含答案) 第一章 概述(包括凸规划) 判断与填空题 arg max f(x)二 arg min 以儿 “ max(x): x D 二 R n 』=-min(x): x D 二 R n ； 设f : D 5 R n > R.若x ： R n ，对于一切R n 恒有f(x”)^f(x)，则称x”为 设f : D 5 R n >R.若x ” ? D ,存在x ”的某邻域N ；(x”)，使得对一切 x ?N .(x)恒有f(x”)：：： f (x)，则称x”为最优化问题 min f (x)的严格局部最 优解? 给定一个最优化问题，那么它的最优值是一个定值 ? V 非空集合D R n 为凸集当且仅当 D 中任意两点连线段上任一点属于 D . V 非空集合D R n 为凸集当且仅当D 中任意有限个点的凸组合仍属于 D . V 任意两个凸集的并集为凸集? 函数f:D R n >R 为凸集D 上的凸函数当且仅当 -f 为D 上的凹函数? V 设f : D R n >R 为凸集D 上的可微凸函数，X ：D ?则对-D ,有 f (x) - f(x )乞 f (x )T (X —X )? 若c(x)是凹函数，则 D={x^R n C(x)启0}是凸集。 V f(x)的算法A 产生的迭代序列，假设算法 A 为下降算法, 则对-k ? 5,1, 2,…匚恒有 ________________ f(x k1)乞 f(x k ) ______________ ? 算法迭代时的终止准则(写出三种) ： ___________________________________________________ 凸规划的全体极小点组成的集合是凸集。 V 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 学号姓名成绩 《冲击式涡轮和反力式涡轮的设计计算》 总结：对冲击式涡轮和反力式涡轮进行设计计算，得到计算结果，具体见表1 和表2。 表1 反力式涡轮的计算结果 表2 冲击式涡轮的计算结果 根据计算结果，我们对比可以得到冲击式涡轮和反力式涡轮的相同点 是： 冲击式涡轮和反力式涡轮在计算功率时，均由泵的功率决定，由 T P N N =∑ 计算。 不同点具体见表3. 表3 反力式涡轮和冲击式涡轮的比较 1. 冲击式涡轮出口压力值取决于涡轮排气是直接排入周围环境还是导入辅助喷管，但两种情况下出口压强和反力式相比均很小。而反力式涡轮通常用于补燃式的液体火箭发动机中的涡轮泵中，所 以在不记喷注器压降的条件下，涡轮的出口压力等于燃烧室的压力。 2.在计算反力式涡轮的参数时，由于反力度容易确定，在分析过程 中广泛采用热力反力度。 反力式涡轮的设计计算 一.反力式涡轮参数的选择 在具有冲击式涡轮的供应系统(无补燃发动机系统)中，由燃气发生器产生的富燃燃气驱动涡轮，涡轮不冷却，富燃燃气的温度在1000～1200K 的范围内，比富氧燃气的允许温度(600～800K)高得多。另外，富燃燃气的气体常数比富氧燃气的气体常数大一些，这些都有利于减小需通过涡轮的燃气流量。 涡轮流量m t q 是具有冲击式涡轮的供应系统的主要参数之一。m t q 值越小，发动机的比冲就越高。涡轮流量m t q 可由泵和涡轮的功率平衡： T Pf Po N N N =+ 泵的需用功率降低，可减小通过涡轮的燃气流量，因此应尽量提高泵的效率。选定泵的结构并确定其效率后，可根据功率平衡求出所需的涡轮燃气流量，由此确定涡轮的效率。 涡轮入口压力(燃气发生器压力)取决于氧化剂泵的出口压力。当用燃料冷却推力室时，燃料泵出口压力比氧化剂泵的出口压力高。 涡轮出口压力之值取决于涡轮排气是直接排入周围环境还是导入辅助喷管。 冲击式涡轮计算的原始数据为： (1)涡轮的设计功率：涡轮功率T N 由泵所需的功率决定，由涡轮泵装置设计任务给定： 其中，T N —涡轮的设计功率，又称涡轮的轴功率； Pf N —燃料泵的轴功率； Pf N —氧化剂泵的轴功率。 (2)涡轮的设计角速度：涡轮的设计转速ω由泵不发生汽蚀时允许的最大角速度确定； (3)涡轮工质的物理常数和温度：涡轮进口总压*0P 、进口总温*0T 、和出口静压2P ；涡轮工质的绝热指数k 和气体常数R 。 二.反力式涡轮参数的选择 一、 填空题 1 ． 若 ()()??? ? ??+???? ?????? ??=212121 312112)(x x x x x x x f ，则 =?)(x f ，=?)(2x f . ２．设f 连续可微且0)(≠?x f ，若向量d 满足 ，则它是f 在x 处的一个下降方向。 3.向量T ) 3,2,1(关于3阶单位方阵的所有线性无关的共轭向量 有 ． 4. 设R R f n →:二次可微，则f 在x 处的牛顿方向为 . 5.举出一个具有二次终止性的无约束二次规划算 法: . 6.以下约束优化问题： )(01)(..)(min 212121 ≥-==+-==x x x g x x x h t s x x f 的K-K-T 条件为： . 7．以下约束优化问题： 1 ..)(min 212 2 21=++=x x t s x x x f 的外点罚函数为(取罚参数为μ) . 二、证明题（７分+8分) 1.设1,2,1,:m i R R g n i =→和m m i R R h n i ,1,:1+=→都是线性函数,证明下 面的约束问题: } ,,1{, 0)(},1{, 0)(..)(min 1112 m m E j x h m I i x g t s x x f j i n k k +=∈==∈≥=∑= 是凸规划问题。 2.设R R f →2 :连续可微，n i R a ∈,R h i ∈,m i ,2,1=,考察如下的约束条件问题: } ,1{,0} 2,1{,0..) (min 11m m E i b x a m I i b x a t s x f i T i i T i +=∈=-=∈≥- 设d 是问题 1 ||||,0,0..)(min ≤∈=∈≥?d E i d a I i d a t s d x f T i T i T 的解，求证：d 是f 在x 处的一个可行方向。 三、计算题（每小题１2分) １.取初始点T x )1,1() 0(=.采用精确线性搜索的最速下降法求解下面的无约束优化问题 (迭代２步）: 2 2212)(m in x x x f += 2.采用精确搜索的BFGS 算法求解下面的无约束问题： 212 2212 1)(min x x x x x f -+= 3.用有效集法求解下面的二次规划问题: . 0,001..42)(min 21212 12 221≥≥≥+----+=x x x x t s x x x x x f 4．用可行方向算法(Zoutend ｉj ｋ算法或Ｆrank Wol ｆｅ算法）求解下面的问题（初值设为)0,0() 0(=x ,计算到)2(x 即可）： . 0,033..22 1)(min 212112 22121≥≥≤+-+-= x x x x t s x x x x x x f 第九章 最优化方法 本章主要介绍线性规划、0-1规划、非线性规划等问题的MATLAB 求解。 9.1 线性规划（Linear Programming ，简写为LP ）问题 线性规划问题就是求多变量线性函数在线性约束条件下的最优值。满足约束条件的解称为可行解，所有可行解构成的集合称为可行域，满足目标式的可行解称为最优解。 MATLAB 解决的线性规划问题的标准形式为： min z f x ￠ =? .. A x b s t Aeq x beq lb x ub ì祝??? ?í??＃??? 其中,,,,,f x b beq lb ub 为列向量，,A Aeq 为矩阵。 其它形式的线性规划问题都可经过适当变换化为此标准形式。 在MATLAB 中求解线性规划问题函数为linprog ，其使用格式为： [x, fval, exitflag, output, lambda] = linprog(f, A, b, Aeq, beq, lb, ub) 输入部分：其中各符号对应线性规划问题标准形式中的向量和矩阵，如果约束条件中有缺少，则其相应位置用空矩阵[]代替。 输出部分：其中x 为最优解，用列向量表示；fval 为最优值；exitflag 为退出标志，若exitflag=1表示函数有最优解，若exitflag=0表示超过设定的迭代最大次数，若exitflag=-2，表示约束区域不可行，若exitflag=-3，表示问题无解，若exitflag=-4，表示执行迭代算法时遇到NaN ，若exitflag=-5，表示原问题和对偶问题均不可行，若exitflag=-7，表示搜索方向太小，不能继续前进；output 表明算法和迭代情况；lambda 表示存储情况。 例1 用linprog 函数求下面的线性规划问题 北京航空航天大学 学年 第一学期期末 《机械设计A4》 考试 A 卷 班 级______________学 号 _________姓 名______________成 绩 _________ 年月日 班号学号姓名成绩 《机械设计A4》考试卷 注意事项： 1、所有题目按步给分，非标准合理答案适当给分，但不超过该步骤的二分之一，计算过程纯计算错误不重复扣分。 2、本试卷共8页，所有题目均在本试题册上作答，拆页或少页本试题册无效。 题目： 一、填空 ……………………………………………………………( 25 分) 二、选择填空 …………………………………………………………( 5 分) 三、简答 ……………………………………………………………( 20 分) 四、分析计算 ……………………………………………………………( 35 分) 五、结构设计 ……………………………………………………………( 15 分) 题号 1 2 3 4 5 成绩 一．填空 ………………………………………………… (共25分，每空0.5分) 1．轴上零件的固定主要是将轴与轴上零件在，和方向上以适当的方式固定。 2．按轴负担的载荷分类，自行车的中轴属于轴；前轴属于轴；后轴 属于轴。 3．带传动的主要失效形式为和，其传动比不稳定主要 是由引起的。 4．闭式软齿面齿轮设计时，考虑到其主要失效形式为 所以一般按 照 强度进行设计，按照 强度进行校核。 5．当滚动轴承在基本额定动载荷作用下运行时，其所能达到的基本额定寿命为 ， 此时滚动轴承的工作可靠度R为。 6．齿轮强度计算中的齿形系数主要取决于 和 。 7．设计中提高轴的强度可以采用、等方法，提高 轴的刚度可以采用等方法。 8．斜齿轮传动与直齿轮相比较，其优点为 、 和 ，开式齿轮传动与闭式齿轮传动比较，其不足之处有 。9．形成流体动力润滑的条件是，， 及。10．三角形螺纹的牙型角α= ，适用于 是因为其 ；矩形 螺纹的牙型角α＝ ，适用于 是因为其 。 11．螺纹防松是要防止 之间的相对运动；常用方法有如，如，如。 12．斜齿轮传动的标准模数是，圆锥齿轮传动的标准模数是， 加工标准直齿轮不发生根切的最小齿数是。 13．代号为71208的滚动轴承，该轴承的类型为,轴承的宽度系列 为，内径尺寸为 mm，精度等级为级。 14．普通平键连接的工作面为，用于轴与轴上零件的固定，传 递。 《结构优化设计》 大作业报告 实验名称: 拓扑优化计算与分析 1、引言 大型的复杂结构诸如飞机、汽车中的复杂部件及桥梁等大型工程的设计问题，依靠传统的经验和模拟实验的优化设计方法已难以胜任，拓扑优化方法成为解决该问题的关键手段。近年来拓扑优化的研究的热点集中在其工程应用上，如: 用拓扑优化方法进行微型柔性机构的设计，车门设计，飞机加强框设计，机翼前缘肋设计，卫星结构设计等。在其具体的操作实现上有两种方法，一是采用计算机语言编程计算，该方法的优点是能最大限度的控制优化过程，改善优化过程中出现的诸如棋盘格现象等数值不稳定现象，得到较理想的优化结果，其缺点是计算规模过于庞大，计算效率太低;二是借助于商用有限元软件平台。本文基于matlab软件编程研究了不同边界条件平面薄板结构的在各种受力情况下拓扑优化，给出了几种典型结构的算例，并探讨了在实际优化中优化效果随各参数的变化，有助于初学者初涉拓扑优化的读者对拓扑优化有个基础的认识。 2、拓扑优化研究现状 结构拓扑优化是近20年来从结构优化研究中派生出来的新分支，它在计算结构力学中已经被认为是最富挑战性的一类研究工作。目前有关结构拓扑优化的工程应用研究还很不成熟，在国外处在发展的初期，尤其在国内尚属于起步阶段。1904 年Michell在桁架理论中首次提出了拓扑优化的概念。自1964 年Dorn等人提出基结构法，将数值方法引入拓扑优化领域，拓扑优化研究开始活跃。20 世纪80 年代初，程耿东和N. Olhoff在弹性板的最优厚度分布研究中首次将最优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的研究兴趣。1988年Bendsoe和Kikuchi发表的基于均匀化理论的结构拓扑优化设计，开创了连续体结构拓扑优化设计研究的新局面。1993年Xie.Y.M和Steven.G.P 提出了渐进结构优化法。1999年Bendsoe和Sigmund证实了变密度法物理意义的存在性。2002 年罗鹰等提出三角网格进化法，该方法在优化过程中实现了退化和进化的统一，提高了优化效率。目前常使用的拓扑优化设计方法可以分为两大类：退化法和进化法。结构拓扑优化设计研究，已被广泛应用于建筑、航天航空、机械、海洋工程、生物医学及船舶制造等领域。 3、拓扑优化建模（SIMP） 结构拓扑优化目前的主要研究对象是连续体结构。优化的基本方法是将设计区域划分为有限单元，依据一定的算法删除部分区域，形成带孔的连续体，实现连续体的拓扑优化。连续体结构拓扑优化方法目前比较成熟的是均匀化方法、变密度方法和渐进结构优化方法。 变密度法以连续变量的密度函数形式显式地表达单元相对密度与材料弹性模量之间的对应关系，这种方法基于各向同性材料，不需要引入微结构和附加的均匀化过程，它以每个单元的相对密度作为设计变量，人为假定相对密度和材料弹性模量之间的某种对应关系，程序实现简单，计算效率高。变密度法中常用的插值模型主要有:固体各向同性惩罚微结构模型(solidisotropic microstructures with penalization，简称SIMP)和材料属性的合理近似模型(rational approximation ofmaterial properties，简称RAMP)。而本文所用即为SIMP插值模型。 1流量工程问题 1.1问题重述 定义一个有向网络G=(N,E)，其中N是节点集，E是弧集。令A是网络G的点弧关联矩阵，即N×E阶矩阵，且第l列与弧里(I,j)对应，仅第i行元素为1，第j行元素为-1， 其余元素为0。再令b m =(b m1 ,…,b mN )T，f m =(f m1 ,…,f mE )T，则可将等式约束表示成： Af m=b m 本算例为一经典TE算例。算例网络有7个节点和13条弧，每条弧的容量是5个单位。此外有四个需求量均为4个单位的源一目的对，具体的源节点、目的节点信息如图所示。这里为了简单，省区了未用到的弧。此外，弧上的数字表示弧的编号。此时，c=((5,5 (5) 1×13 )T， 根据上述四个约束条件，分别求得四个情况下的最优决策变量x=((x 12,x 13 ,…,x 75 ) 1×13 )。 图 1 网络拓扑和流量需求 1.27节点算例求解 1.2.1\ T） 1.2.2算例1（b1=[4;-4;0;0;0;0;0] 转化为线性规划问题： Minimize c T x1 Subject to Ax1=b1 x1>=0利用Matlab编写对偶单纯形法程序，可求得: 最优解为x1*=[4 0 0 0 0 0 0 0 0 0 0 0 0]T 对应的最优值c T x1=20 1.2.3算例2（b2=[4;0;-4;0;0;0;0]T） Minimize c T x2 Subject to Ax2=b2 \ X2>=0利用Matlab编写对偶单纯形法程序，可求得: 最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T 对应的最优值c T x2=20 1.2.4算例3（b3=[0;-4;4;0;0;0;0]T） Minimize c T x3 Subject to Ax3=b3 X3>=0利用Matlab编写对偶单纯形法程序，可求得: 最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T最优化方法简明教程—centre
北航最优化方法大作业参考
最优化方法大作业答案
最优化方法及应用
最优化方法大作业
北航计算机控制系统大作业
北航惯性导航大作业
北航航空工程大型通用软件应用大作业样本
北航计算机控制系统大作业
最优化计算方法课后习题答案----高等教育出版社。施光燕
北航数值分析大作业第二题精解
《最优化方法》复习题(含答案)
北航涡轮泵大作业
最优化方法(试题+答案)
第九章 最优化方法
北航机械设计试题
结构优化设计大作业(北航)
北航最优化方法大作业参考