北方工业大学建模与仿真

1、系统建模与仿真2011.11.151第一节：基本概念实际系统计算机仿真系统数学模型2目的3目的4目的5目的，应用，操作训练数字计算试验军事领域: 武器装备,军事训练,先进概念工业领域:核电站,电力工业,虚拟制造技术教育训练: 载体操作,过程控制,博弈决策民用工程:交通流仿真,输水工程仿真,人口仿真,疾病扩善仿真等等 6回顾电子计算机的出现应用于航空航天领域工业操作仿真社会学仿真半物理仿真集成化建模与仿真环境训练,视觉,虚拟等等仿真技术 7发展趋势 面向部件的系统仿真工具智能仿真：智能，人脑活动仿真，行为仿真可视化仿真技术：3维视觉虚拟现实仿真：数据手套，数据眼镜，数据头盔等等，用于科学试验。I

2、nternet仿真：新一代的仿真试验，在Internet上试验，多人参与的、远程的、交互的试验。8本课程讲解内容连续系统建模与数值仿真 前8节四旋翼机器人建模及程序设计 后8节9举例1：人体运动仿真 10举例2：电池组冷却气流仿真 11举例3：航空发动机仿真12仿真环境13输出曲线14四旋翼机器人-视频15自适应动态逆控制16第二节：连续系统的微分方程建模x 矢量t 时间f 函数动态系统17数值求解1：显式欧拉方法 k 上标k代表第k步 时间为t=khh 代表时间步长 k到k+1步的时间间隔 显示格式：直接求解时间步长：受约束计算精度：1阶精度18计算稳定性19数值求解2：隐式欧拉方法 k 上

3、标k代表第k步 时间为t=khh 代表时间步长 k到k+1步的时间间隔 显示格式：求解非线性方程租时间步长：不受约束计算精度：1阶精度20计算稳定性21 改进欧拉方法精度问题 二阶精度，计算精度；显式方法，受时间步长制约22 四阶龙格-库塔法精度问题 简略写出，显式格式，受时间步长制约23微分方程数值求解代码： 欧拉方法： (*fty)(y,d,n); for(i=0;i=n-1;i+)yi=yi+h*di; 改进欧拉方法： for(i=0;i=n-1;i+)zi=yi; (*fty)(y,d,n); for(i=0;i=n-1;i+)yi=zi+h*di; (*fty)(y,d,n); fo

4、r(i=0;i=n-1;i+)di=zi+h*di; for(i=0;i=n-1;i+)yi=(yi+di)/2;龙哥库塔方法：for(i=0;i=n-1;i+)zi=yi;bi=yi; (*fty)(y,d,n); for(j=0;j=2;j+) for(i=0;i=n-1;i+) yi=zi+aj*di; bi=bi+aj+1*di/3.0; /s=t+aj; (*fty)(y,d,n); for(i=0;iU) Uc=U; endif(Uceps G(fk),J(fk)fk+1=f k +f fkCriterion f is okf*60Case A: component fault 6

5、1第六节 神经网络建模有些时候我们对一个系统的工作机理不是很清楚的时候，我们无法获得该系统的解析模型，就是无法得到系统的微分方程模型。这时候，我们往往能够得到系统的一些输入输出关系数据，那么根据输入输出关系数据建立的系统模型，我们成为黑箱模型。建立黑箱模型中最为经典模型为神经网络模型。建模过程就是在寻找一种输入变量序列与输出变量序列的一种数学映射关系，期望找到一个通用的模型，神经网络就是其中的一种通用模型。目前已经研制出包括bp,rbf,fir,art, Kohonen等数百种网络以及各种经典算法。起源于对生物神经元的研究，生物神经元的基本构成：突触，核，树突等，互相连接。径向基神经网络建立系

6、统模型，这种网络的优点是泛化能力强，不需要复杂的学习过程，工程应用价值较高。 62径向基函数神经网络 63数学表达与模型的由来第一层计算模型： 第二层计算模型宽度为一维参数，小非线性程度，大去除数据噪声广义非线性回归过程，无限拟合观测数据64 学习算法之一：最近邻聚类选择一个适当的高斯函数的宽度从第一个数据序列开始，建立的RBF网络只有一个隐含单元 考虑第k个样本序列，存在M个聚类点，求出相对于M个聚类中心点的距离 聚类原则： 小于高斯函数宽度聚类，否则新建聚类中心65应用举例：股指预测股票价格预测便受到广泛关注与积极研究。股票价格是多变量的动态非线性系统，影响因素众多，包括历史数据、技术指标

7、、宏观因素等等各种变量，只有从中找出一组适当的影响因素，才能有效地解释股票市场价格的变化关系 技术分析预测价格走势中，K线数据是所有技术指标的基础，K线数据包括：开盘价、最高价、最低价、收盘价和成交量，5个数据构成时间序列，那么就可以根据最近几个交易日的K线数据预测下一个交易日的K线数据。在对K线数据进行无量纲化处理后，就可以采用非线性预测模型预测股票价格的未来走势。66股指数据结构及导数处理 收盘价开盘价最高价最低价成交量收盘价导数=（收盘价-前一日收盘价）/前一日收盘价最高价导数=（最高价-前一日收盘价）/前一日收盘价最低价导数=（最低价-前一日收盘价）/前一日收盘价开盘价导数=（开盘价-

8、前一日收盘价）/前一日收盘价成交量导数=（成交量-前一日成交量）/前一日成交量 67K线形态数据的预先处理 处理前处理后如果选取第i个交易日的前w个交易日的K线形态数据为模型输入，选取第i+1个K线形态数据为模型输出，则有约束为i+1小于等于n, i大于w,则称w为预测窗口宽度 68径向基神经网络的泛化学习：聚类其中为 隐层单元第l个结点高斯函数中心，由聚类方法产生69沪深300指数实例预测 日前K线无量纲数据和其后预测的67个交易日无量纲数据 日前成交量无量纲数据和其后预测的67个交易日无量纲数据 70沪深300指数实例预测日前K线形态和其后预测的67个交易日形态 日前成交量形态和其后预测的

9、67个交易日形态 71建模与应用72所谓指数乖离法73神经网络预测举例74建模与仿真演讲作业1如上图所示直流电机电路图，电机转动惯量为J,电机扭矩方程为，转子电流方程为：，定子电流方程，回答：将如上方程写为的状态空间模型形式。求解在Ua= ？ ,Ue= ？ 输入条件下，系统稳定的转速和ia,ie大小？在2稳定条件下，若Ua在？秒内变为？，Ue在？秒内变为？，试求10s内的转子转速的动态数值解，并给出转子电流和定子电流数值解？75建模与仿真演讲作业21.假设系统可观测量为位置x和速度u，将车辆跟驰模型写为关于位置x,速度u的状态空间模型2.写出你自己的车辆自由方程u(xn-xn+1)3.假设测量

10、得到n车的速度和位置，即xn和un, ,编写后车的系统响应数值解76总结状态方程：非线性，状态量求解状态方程方法：原生态数值解与解析解对向量的理解问题77四旋翼机器人建模78多轴旋翼机器人四旋翼机器人-Quadrotor Robot多轴-四轴基础上多轴四轴:1对正转-1对反转 十字排列动力:无刷电机带动螺旋浆动作:垂直起飞-着陆-前进-后退-侧飞-转身-倒飞特点:飞行更灵活,任意方向飞行,噪声小,狭小空间起飞着陆,适应性和隐蔽性更强缺点:受电池供电限制，其航时短，有效载荷小，任务受限制，其控制系统复杂79应用领域广告:广告宣传，其小巧，携带方便，飞行成本低廉，独特造型吸引客户群，可以完成空中拍

11、照等任务娱乐:航模比赛，吸引更多的航模爱好者，航空摄影，为影视公司提供可操控性和稳定性服务工业:电力巡线、电力线路检测、电力架线，由于其体积小，使用方便，可以进入受设备大小无法涉足的领域完成巡查或则检测任务军事:适用战地侦查、反恐，其噪音小、隐蔽性强，未来可以广泛应用在山地作战、城市作战反恐等领域80多轴旋翼机器人81国内研究现状高校院所:北京航空航天大学，南京航空航天大学，南京理工大学，哈尔滨工业大学，中国科学技术大学，东北大学等高校先后投入人力和物力，进行算法、控制、工程方面的研究科技公司:国内一些科技公司也相继开发多种民用型四旋翼或多旋翼飞行器，用于监控、侦查、航空摄影等名种任务，目前开

12、发和研究刚刚起步阶段,如XAirCraft等航模爱好者:由于四旋翼机器人其结构简单，机动性好，控制灵活，未来有很大的军事和民用发展空间吸引大批航模爱好者 82刚体旋翼机器人建模平地坐标系 NED North-East-Down 站立向北看右手坐标系,习惯上为正方向,欧拉角:机体坐标系Body Co-Ordinate 右手坐标系,体轴坐标原点与质心一致姿态:体坐标系与平地坐标系之间的关系欧拉角: 滚转角:0-180度 俯仰角:0-90度 方位角:0-360度,正北:0度 东:90度 zN83刚体动力学-Newton-Euler方程Mass of rigid bodySum of TorqueSu

13、m of Force(1)Force equation (2)Moment equation intertia matrix84导航(3)Rotation transform(4)Navigation transformEuler anglePosition 85外力及外力矩分析Distance between motor and coordinate center LAeroforce and aerotorque coefficientX type configuration86螺旋浆拉力、陀螺力模型i=1，2，3，4螺旋浆角速度拉力转换系数力矩转换系数简化： 1) 刚体运动，质量不变 2

14、) 平地假设，重力加速度不随高度改变 3) 机体形状和质量对称 *4) 速度低，近似气动力与速度有关87无刷直流电机转动力学模型电枢电流，供电电压，电机电感，电机内阻，反电势系数电机角速度，电机及负载转动惯量实际建模时，考虑四旋翼无刷电机的电感很小，可忽略电流动力学方程摩擦系数，电机转矩系数，载荷系数88旋翼机器人的动力学模型(1)Force equation (2)Moment equation(3)Rotation transform(4)Navigation transform(5)Electric motor dynamic89动力学模型状态量:12+4=16个控制量:4个控制电压参数

15、集:6个拉力转换系数力矩转换系数悬臂梁柔性系数摩擦系数，电机转矩系数，载荷系数气动力和力矩系数90四旋翼飞行器控制与仿真91电机动力学parameterdescriptionValue in quadrotor dynamic modelmodeling error percent in adaptive dynamic inversionmMass of quadrotor2.3kg120%Lr,Lx,LyLever arm length0.35m,0.247m,0.247m100%IxxBody inertia of x axis8.04*10-3 kgm275%IyyBody inert

16、ia of y axis8.46*10-3 kgm2120%IzzBody inertia of z axis14.68*10-3 kgm280%kfForce coefficient0.65016e-3N/(rad/s)295%ktTorque coefficient0.82218e-5N/(rad/s)285%k1,k2,k3Motor model coefficient20, 0.01, 3.50klElastic coefficient of gear6000N/m092四旋翼机姿态与轨迹的线性方法时间尺度分析：电机控制（极快回路），快回路（姿态控制），慢回路（轨迹控制）93PD控制器

17、,94自适应动态逆控制自适应项95Lyapunov稳定性的渐进稳定性96仿真程序设计function dxv=quad(xv,uv,kv)Mass=2.3 ; %kgLeverX=0.60 ; %meterLeverY=0.35 ; %meterIxx=8.04*10-3 ; %kgm2Iyy=8.46*10-3 ; %kgm2Izz=14.68*10-3 ; %kgm2kl=6000 ;zl=Mass*9.8/kl ;Ir=5.861*10-5 ; %propeller gyroscoptic inertia ub=xv(1) ; %m/svb=xv(2) ; %m/swb=xv(3) ;

18、%m/spb=xv(4) ; %rad/sqb=xv(5) ; %rad/srb=xv(6) ; %rad/sphi=xv(7) ; %radthe=xv(8) ; %radpsi=xv(9) ; %radxe=xv(10) ; %mye=xv(11) ; %mze=xv(12) ; %m %define the state varibale of four motorw1=xv(13) ; %rad/sw2=xv(14) ; %rad/sw3=xv(15) ; %rad/sw4=xv(16) ; %rad/s u1=uv(1) ; %u1 1th motoru2=uv(2) ; %u2 2t

19、h motoru3=uv(3) ; %u3 3th motoru4=uv(4) ; %u4 4th motor k1=kv(1) ;k2=kv(2) ;k3=kv(3) ;97%define the propel force kf=kv(4) ;kt=kv(5) ;Mp=kf*LeverX*(-w1*w1 + w2*w2 + w3*w3 - w4*w4) ;Mq=kf*LeverY*(-w1*w1 + w2*w2 - w3*w3 + w4*w4) ;Mr=kt*(LeverX2+LeverY2)0.5*( w1*w1 + w2*w2 - w3*w3 - w4*w4) ;Fz=kf* ( w1*

20、w1 + w2*w2 + w3*w3 + w4*w4) ; %define the ground support or elastic force of body landing if (zl-ze)0 Fl=kl*(zl-ze) ;else Fl=0 ;end%define body motion kub=0.2;kvb=0.2;kwb=0.25;dub= -( -rb*vb + qb*wb) + 9.8*sin(the) + kub*ub/Mass;dvb= -( rb*ub - pb*wb) - 9.8*cos(the)*sin(phi) + kvb*vb/Mass;dwb= -(-qb

21、*ub +pb*vb ) - 9.8*cos(the)*cos(phi) + Fz/Mass + Fl/Mass + kwb*wb/Mass; %dpb= - ( -rb*Iyy*qb +qb*Izz*rb)/Ixx + Mp/Ixx%dqb= - ( rb*Ixx*pb -pb*Izz*rb)/Iyy + Mq/Iyy%drb= - (-qb*Ixx*pb +pb*Iyy*qb )/Izz + Mr/Izzkpb=0.00001; kqb=0.00001; krb=0.00001;dpb= - (-Iyy + Izz)/Ixx*qb*rb + Mp/Ixx + Ir*(w1+w2-w3-w4

22、)*qb/Ixx + kpb*pb/Ixx;dqb= - ( Ixx - Izz)/Iyy*pb*rb + Mq/Iyy - Ir*(w1+w2-w3-w4)*pb/Iyy + kqb*qb/Iyy;drb= - (-Ixx + Iyy)/Izz*pb*qb + Mr/Izz + 0 + krb*rb/Izz; 98%define the attitude dynamicdphi= pb + qb*sin(phi)*tan(the) + rb*cos(phi)*tan(the) ; dthe= qb*cos(phi) - rb*sin(phi) ;dpsi= qb*sin(phi)/cos(t

23、he) + rb*cos(phi)/cos(the) ; %define navigation dynamicdxe= ub*cos(the)*cos(psi) + vb*(-cos(phi)*sin(psi)+sin(phi)*sin(the)*cos(psi) + wb*( sin(phi)*sin(psi)+cos(phi)*sin(the)*cos(psi);dye= ub*cos(the)*sin(psi) + vb*(+cos(phi)*cos(psi)+sin(phi)*sin(the)*sin(psi) + wb*(-sin(phi)*cos(psi)+cos(phi)*sin

24、(the)*sin(psi);dze= ub*(-sin(the) + vb*sin(phi)*cos(the) + wb*cos(phi)*cos(the) ; dw1= -k1*w1 -k2*w1*w1 +k3*u1 ;dw2= -k1*w2 -k2*w2*w2 +k3*u2 ;dw3= -k1*w3 -k2*w3*w3 +k3*u3 ;dw4= -k1*w4 -k2*w4*w4 +k3*u4 ; %derivertive of state vectordxv=dub dvb dwb dpb dqb drb dphi dthe dpsi dxe dye dze dw1 dw2 dw3 dw

25、499慢回路控制function uv=slowloop(xv,xcmd,deltamass)xe=xv(1) ; %mye=xv(2) ; %mze=xv(3) ; %mpsi=xv(4) ; %rad xe1=xv(5) ;%approximation x axis velocity in earth frameye1=xv(6) ;%approximation y axis velocity in earth frameze1=xv(7) ;%approximation z axis velocity in earth frame %define the refference comma

26、ndxecmd =xcmd(1)yecmd =xcmd(2)zecmd =xcmd(3)psicmd=xcmd(4) Mass=2.3*1.2 ; %kg %define the outer loop of inverse dynamicksi=0.72 %damping ratioomiga=1.0 %natural errxe=xe-xecmd; errxe=positionlimit(errxe); errye=ye-yecmd; errye=positionlimit(errye);errze=ze-zecmd; errze=positionlimit(errze); xe2 = -2

27、*ksi*omiga*xe1 - omiga*omiga*errxe ;ye2 = -2*ksi*omiga*ye1 - omiga*omiga*errye ;ze2 = -2*ksi*omiga*ze1 - omiga*omiga*errze ; deltamass = deltamass + 0.2*ze2; Fz=(ze2 + 0.7*deltamass + 9.8 )*Mass ; %deltamass :adaptive sliding mode dynamic inversion controlphicmd=(xe2*sin(psi)-ye2*cos(psi)*Mass/Fz ;%

28、45 degree limit for roll angleif phicmd0.7854 phicmd=0.7854;end%45 degree limit for pitch anglethecmd=(xe2*cos(psi)+ye2*sin(psi)*Mass/Fz ;if thecmd0.7854 thecmd=0.7854;end psicmd=psicmd; uv=phicmd thecmd psicmd Fz deltamass; %command for fastloop100快回路控制function uv=fastloop(xv,xcmd,deltaangle)phi=xv

29、(1) ; %radthe=xv(2) ; %radpsi=xv(3) ; %radpb=xv(4) ; %rad/sqb=xv(5) ; %rad/srb=xv(6) ; %rad/s%define the refference commandphicmd =xcmd(1)thecmd =xcmd(2)psicmd =xcmd(3)deltapb=deltaangle(1)deltaqb=deltaangle(2)deltarb=deltaangle(3) errphi=(phicmd-phi); errphi=anglelimit(errphi);errthe=(thecmd-the);

30、errthe=anglelimit(errthe);errpsi=(psicmd-psi); errpsi=anglelimit(errpsi);%define the parameter of modelIxx=8.04*10-3*0.75 ; %kgm2Iyy=8.46*10-3*1.2 ; %kgm2Izz=14.68*10-3*0.8 ; %kgm2kphi=7.14 ; %5kthe=7.14 ; %5kpsi=1.43 ; %5R3=1 sin(phi)*tan(the) cos(phi)*tan(the);0 cos(phi) -sin(phi); 0 sin(phi)/cos(

31、the) cos(phi)/cos(the)v3=kphi*errphi kthe*errthe kpsi*errpsiv3=inv(R3)*v3pbcmd=v3(1)qbcmd=v3(2)rbcmd=v3(3) %define the inner loop for attitude augment systemkpb=14;%30;%25kqb=14;%30;%25krb=2.8;%30;%25 errpb = pb - pbcmd;errqb = qb - qbcmd;errrb = rb - rbcmd; %adaptive variable deltaangle dynamic pro

32、cessdeltapb = deltapb + 0.02*errpb;deltaqb = deltaqb + 0.02*errqb;deltarb = deltarb + 0.02*errrb;%v3=-kpb*errpb*Ixx -kqb*errqb*Iyy -krb*errrb)*Izz;v3=-kpb*(errpb+ 0.7*deltapb)*Ixx -kqb*(errqb+0.7*deltaqb)*Iyy -krb*(errrb+0.17*deltarb)*Izz;v3=v3 + (Izz-Iyy)*qb*rb (Ixx-Izz)*pb*rb (Iyy-Ixx)*pb*qb; Tp=v

33、3(1);Tq=v3(2);Tr=v3(3); uv=Tp Tq Tr deltapb deltaqb deltarb;101逆变换function uv=inverse(xv)Tp=xv(1)Tq=xv(2)Tr=xv(3)Fz=xv(4)%define the parameter of modelLever=0.35 ; %meterkf=0.650157e-3*0.95kt=0.82218e-5*0.85 %define the motor dynamic commandR4=-1 1 1 -1; -1 1 -1 1; 1 1 -1 -1; 1 1 1 1v4=inv(R4)*Tp/(k

34、f*Lever) Tq/(kf*Lever) Tr/(kt*Lever) Fz/kf v4(1)=ceilfloor(v4(1)v4(2)=ceilfloor(v4(2)v4(3)=ceilfloor(v4(3)v4(4)=ceilfloor(v4(4) ksi=1000./163.8; %4s Li batteryuv=sqrt(v4)*ksi;102主程序及初始化function inversion() ub=0;vb=0;wb=0;pb=0;qb=0;rb=0;phi=0;the=0;psi=0;xe=0;ye=0;ze=0;w1=0;w2=0;w3=0;w4=0; deltamass=

35、0; %adaptive variable of delta mass gravity deltapb=0; %adaptive variable of delta roll angle phideltaqb=0; %adaptive variable of delta pitch angle psideltarb=0; %adaptive variable of delta yaw angle psi xv=ub vb wb pb qb rb phi the psi xe ye ze w1 w2 w3 w4; %the state for quadrotor kf=0.650157e-3;k

36、t=0.82218e-5;timeresponse=0.1;voltage=25.2;speed=163.8kv=20 0.01 3.5 kf kt; timestep=0navigationpx=0navigationpy=0navigationpz=0direction=0 sw1=0;sw2=0;sw3=0;sw4=0;delta1=0;delta4=0;delta5=0;delta6=0; trace= 0 0 10 0; 10 0 10 0; 10 10 10 0; 0 10 10 0; 0 5 10 0; 0 5 10 3; %trace= 0 0 10 0;% 1 1 10 0.

37、5; xstep=xe;ystep=ye;zstep=ze; 103仿真主程序循环for kstep=1:6 xcmdslow=trace(kstep,:) % if i=25% xcmdslow=1 2 2 0% end for i=1:60 xslow=xe ye ze psi (xe-xstep)/0.2 (ye-ystep)/0.2 (ze-zstep)/0.2 xcmdfast=slowloop(xslow,xcmdslow,deltamass) %adaptive variable mass gravitty xstep=xe; ystep=ye; zstep=ze; Fz=xcm

38、dfast(4) deltamass= xcmdfast(5) %the fastloop for loopvariable k for k=1:10 xfast=phi the psi pb qb rb deltaangle=deltapb deltaqb deltarb; ucmd=fastloop(xfast,xcmdfast,deltaangle) %adaptive variable for atitude angle Tp=ucmd(1) Tq=ucmd(2) Tr=ucmd(3) deltapb=ucmd(4) deltaqb=ucmd(5) deltarb=ucmd(6) uv

39、cmd=Tp Tq Tr Fz uv=inverse(uvcmd) for m=1:10 xv=xv+0.002*quad(xv,uv,kv) end ub=xv(1);vb=xv(2);wb=xv(3); pb=xv(4);qb=xv(5);rb=xv(6); phi=xv(7);the=xv(8);psi=xv(9); xe=xv(10);ye=xv(11);ze=xv(12); w1=xv(13);w2=xv(14);w3=xv(15);w4=xv(16); timestep=timestep (kstep-1)*60*0.2+(i-1)*0.2+k*0.02; navigationpx

40、=navigationpx xe; navigationpy=navigationpy ye; navigationpz=navigationpz ze; direction=direction psi; sw1=sw1 phi;sw2=sw2 the;sw3=sw3 psi;sw4=sw4 ze; delta1=delta1 deltamass; delta4=delta4 deltapb;delta5=delta5 deltaqb;delta6=delta6 deltarb; end end end 104姿态与轨迹绘图部分 figure plot(timestep,navigationp

41、x,timestep,navigationpy,timestep,navigationpz,timestep,direction) title(position coordinate trajectory tracking output) figure plot(timestep,sw1,timestep,sw2,timestep,sw3,timestep,sw4) title(attitude angle output) figure plot(timestep,delta1,timestep,delta4,timestep,delta5,timestep,delta6) title(ada

42、ptive variable dynamic output) figure plot3(navigationpx,navigationpy,navigationpz) title(trajectory seeking)105最优化计算方法 无约束最优化问题的定义 min f(x) s.t. x X多元函数f(x) 106二次函数的梯度和Hesse矩阵107无约束最优化条件存在0，使满足所有x X 且|x-x*|=f(x*),则称x*为f的局部极小点若f(x)连续可微，x*为局部极小点，则g(x*)=0若f(x)二阶连续可微，x*为局部极小点的充分条件g(x*)=0,H(x*)是正定或半正定矩阵

43、108无约束优化问题的算法结构第k步位移第k步步长第k步搜索方向基本步骤：a)给定初始点x0 b)确定xk的下降方向dk c)确定步长k,使f(xk+ kdk)0,k=0 b)检验终止条件，计算dk=-gk,若| gk |0,x0,k=0 b)检验终止条件，计算gk,| gk | err x=x-alpha*g g=dfun(x) d=sqrt(g(1)2+g(2)2)end function g=dfun(x)h=0.01j=fun(x)for i=1:2 z=x(i) x(i)=x(i) + h dj=fun(x) g(i)=(dj-j)/h x(i)=zEndfunction J=fun

44、(x)J=4*x(1)2+x(2)2-x(1)2*x(2)113err=0.0001alpha=0.00002x=-1.2,1 g=dfun1(x)d=sqrt(g(1)2+g(2)2)while d err x=x - alpha*g g=dfun1(x) d=sqrt(g(1)2+g(2)2)end function g=dfun1(x)h=0.00001;j=fun1(x);for i=1:2 z=x(i); x(i)=x(i) + h; dj=fun1(x); g(i)=(dj-j)/h; x(i)=z;Endfunction J=fun1(x)J=100*(x(2)-x(1)2)2

45、+ (1-x(1)2;114拟牛顿法最有效的无约束优化方法模拟Hesse矩阵，根据1阶导数信息及位移，构造二阶导数矩阵相似正定矩阵拟牛顿条件115拟牛顿法构造Hesse矩阵逆116DFP方法计算步骤：a)选定初值x0,误差e,H0对称正定，k=0 b)检验终止条件，计算gk,若|gk| err sk=-alpha*inv(Hk)*g x=x+sk g1=dfun(x) yk=g1-g vk=sqrt(yk*Hk*yk)*(sk/(sk*yk)-Hk*yk/(yk*Hk*yk) Hk=Hk + sk*sk/(sk*yk) - Hk*yk*yk*Hk/(yk*Hk*yk) + vk*vk g=g1

46、 d=sqrt(g(1)2+g(2)2)end 119信赖域方法优化步长位移，直接确定步长位移首先给出位移长度上界，以当前点为中心，以上界为半径确定领域，然后通过求该领域的二次近似目标函数的最优点来确定候选位移，此区域内近似函数往往是可信赖的而称为信赖域方法。信赖域半径需要调整以使目标函数充分下降120实现步骤选初值x0,选信赖域半径0,误差e检验终止条件，|gk|，有x k+1 =x k +s k,否则x k+1 =x k121如何求解子问题Cauchy MethodDogleg Method122信赖域法的计算程序err=0.0001et1=0.01et2=0.75gm1=0.5gm2=2.0 x=1;1 g=dfun(x)d=sqrt(g(1)2+g(2)2)delta=d/10Hk=1 0;0 1while d err tau=-delta/d*g if(tau0) sk=-delta/d*g else sk=-min(d2/(g*Hk*g),delta/d)*g endfk=fun(x) Ar=fk-fun(x+sk) Pr=fk+g*sk+sk*Hk*sk/2 Rk=Ar/Prif (Rk=et2) delta=gm2*delta else delta=delta end