最优化方法大作业-算法源程序-0.618法、抛物线法、共轭梯度法

最优化方法大作业-算法源程序-0.618法、抛物线法、共轭梯度法
最优化方法大作业-算法源程序-0.618法、抛物线法、共轭梯度法

最优化方法程序作业

专业:油气储运工程

班级:

姓名:

学号:

一、开发工具

该应用程序采用的开发工具是Visual studio 2005,编程语言使用的是C#。

二、程序代码

(1)0.618法和抛物线法求ψ(t)=sint,初始点t0=1

①主程序:

using System;

using System.Data;

using System.Configuration;

using System.Web;

using System.Web.Security;

using System.Web.UI;

using System.Web.UI.WebControls;

using System.Web.UI.WebControls.WebParts;

using System.Web.UI.HtmlControls;

using System.Collections;

public partial class_Default : System.Web.UI.Page

{

JiSuan JS = new JiSuan();

protected void Button1_Click(object sender, EventArgs e)

{

double xx0 = 0.01;//步长

double tt0 = 1, pp = 2, qq = 0.5;

ArrayList list1 = new ArrayList();

list1 = (ArrayList)JS.SuoJian(xx0, tt0, pp, qq);//调用倍增半减函数double aa = double.Parse(list1[0].ToString());

double bb = double.Parse(list1[1].ToString());

txtShangxian.Text = bb.ToString();//在页面上显示极小区间

txtXiaxian.Text = aa.ToString();

ArrayList list2 = new ArrayList();

list2 = (ArrayList)JS.JiXiao1(aa, bb);//调用0.618法函数

double jixiao1 = double.Parse(list2[0].ToString());

double fjixiao1 = double.Parse(list2[1].ToString());

txtJixiao1.Text = jixiao1.ToString();//在页面上显示极小点

txtFjixiao1.Text = fjixiao1.ToString();//在页面上显示极小点处的函数值

ArrayList list3 = new ArrayList();

list3 = (ArrayList)JS.JiXiao2(aa, bb);//调用抛物线法函数

double jixiao2 = double.Parse(list3[0].ToString());

double fjixiao2 = double.Parse(list3[1].ToString());

txtJixiao2.Text = jixiao2.ToString();//在页面上显示极小点

txtFjixiao2.Text = fjixiao2.ToString();//在页面上显示极小点处的函数值 }

}

②各个子函数的程序:

using System;

using System.Data;

using System.Configuration;

using System.Web;

using System.Web.Security;

using System.Web.UI;

using System.Web.UI.WebControls;

using System.Web.UI.WebControls.WebParts;

using System.Web.UI.HtmlControls;

using System.Collections;

///

/// JiSuan 的摘要说明

///

public class JiSuan

{

public JiSuan()

{

//

// TODO: 在此处添加构造函数逻辑

//

}

//目标函数

public double f(double z)

{

double zz = Math.Sin(z);

return zz;

}

//倍增半减函数

public ArrayList SuoJian(double x0, double t0, double p, double q)

{

double t1;

double f0, f1;

double x1, t2;

double f2;

double a = 0, b = 0, c = 0;

double fa, fc, fb;

double t3, f3;

bool biaozhi1 = true;//设置是否循环的标志

bool biaozhi2 = true;

t1 = t0 + x0;

f0 = f(t0);//调用目标函数

f1 = f(t1);

if (f0 > f1)

{

while (biaozhi1)

{

x1 = p * x0;

t2 = t1 + x1;

f2 = f(t2);//调用目标函数

if (f1 > f2)

{

t0 = t1;

t1 = t2;

f0 = f1;

f1 = f2;

continue;

}

else

{

a = t0;

c = t1;

b = t2;

fa = f0;

fc = f1;

fb = f2;

break;

}

}

}

else

{

t3 = t1;

f3 = f1;

t1 = t0;

f1 = f0;

t0 = t3;

f0 = f3;

while (biaozhi2)

{

x1 = q * x0;

t2 = t1 - x1;

f2 = f(t2);//调用目标函数

if (f1 > f2)

{

t0 = t1;

t1 = t2;

f0 = f1;

f1 = f2;

continue;

}

else

{

a = t2;

c = t1;

b = t0;

fa = f2;

fc = f1;

fb = f0;

break;

}

}

}

ArrayList list = new ArrayList();

list.Add(a);

list.Add(b);

return list;

}

//0.618法函数

public ArrayList JiXiao1(double a, double b) {

double jd = 0.0001;//精度

double p = 0.618;

double t1, t2;

double f1, f2;

double jixiao=0;//极小点

bool biaozhi1 = true;//设置是否循环标志

bool biaozhi2 = true;

while (biaozhi1)

{

t1 = a + (1 - p)*(b - a);

t2 = a + p * (b - a);

f1 = f(t1);//调用目标函数

f2 = f(t2);

while (biaozhi2)

{

if (f1 < f2)

{

b = t2;

t2 = t1;

f2 = f1;

t1 = a + (1 - p)*(b - a); f1 = f(t1);//调用目标函数

if (Math.Abs(b - a) > jd) {

continue;

}

else

{

biaozhi1 = false;

break;

}

}

else if (f1 == f2)

{

a = t1;

b = t2;

if (Math.Abs(b - a) > jd) {

break;

}

else

{

biaozhi1 = false;

break;

}

}

else if (f1 > f2)

{

a = t1;

t1 = t2;

f1 = f2;

t2 = a + p * (b - a);

f2 = f(t2);//调用目标函数

if (Math.Abs(b - a) > jd)

{

continue;

}

else

{

biaozhi1 = false;

break;

}

}

}

}

jixiao = (a + b) / 2;

double fjixiao = f(jixiao);//调用目标函数

ArrayList list = new ArrayList();

list.Add(jixiao);

list.Add(fjixiao);

return list;

}

//抛物线法函数

public ArrayList JiXiao2(double a, double b)

{

double jd = 0.0001;//精度

double c = a + (b - a) / 2;//将c的值设为a、b点的中点

double fa, fb, fc, c1, c2;

double jixiao = 0;//极小点

double ta, fta;

double fac;

bool biaozhi1 = true;//设置是否循环标志

while (biaozhi1)

{

fa = f(a);//调用目标函数

fb = f(b);

fc = f(c);

c1 = (fb - fa) / (b - a);

c2 = ((fc - fa) / (c - a) - c1) / (c - b);

if (c2 == 0)

{

jixiao = c;

break;

}

else

{

ta = 0.5 * (a + b - (c1 / c2));

fta = f(ta);//调用目标函数

if (Math.Abs(b - a) <= jd)

{

jixiao = ta;

break;

}

else

{

if (fc > fta)

{

if (c > ta)

{

b = c;

fb = fc;

c = ta;

fc = fta;

continue;

}

else

{

a = c;

fa = fc;

c = ta;

fc = fta;

continue;

}

}

else if (fc == fta)

{

if (c > ta)

{

a = ta;

b = c;

c = (c + ta) / 2;

fa = fta;

fb = fc;

fc = f(c);//调用目标函数

continue;

}

else if (c == ta)

{

fac = f((a + c) / 2);//调用目标函数

if (fac < fc)

{

c = (a + c) / 2;

fc = f(c);//调用目标函数

b = ta;

fb = fta;

continue;

}

else

{

a = (a + c) / 2;

fa = f(a);//调用目标函数

continue;

}

}

else if (c < ta)

{

a = c;

c = (c + ta) / 2;

b = ta;

fa = fc;

fb = fta;

fc = f(c);//调用目标函数

continue;

}

}

else if (fc < fta)

{

if (c > ta)

{

a = ta;

fa = fta;

continue;

}

else

{

b = ta;

fb = fta;

continue;

}

}

}

}

}

double fjixiao = f(jixiao);//调用目标函数

ArrayList list = new ArrayList();

list.Add(jixiao);

list.Add(fjixiao);

return list;

}

}

(2)共轭梯度法求函数极小点:f(x)=1.5x12+0.5x22-x1x2-2x1;初始点为(-2,4)①主程序:

using System;

using System.Data;

using System.Configuration;

using System.Web;

using System.Web.Security;

using System.Web.UI;

using System.Web.UI.WebControls;

using System.Web.UI.WebControls.WebParts;

using System.Web.UI.HtmlControls;

using System.Collections;

public partial class_Default : System.Web.UI.Page

{

JiSuan JS = new JiSuan();

protected void Button1_Click(object sender, EventArgs e)

{

//共轭梯度法求极小值

double[] x0 = new double[2];//定义二维数组

double[] x = new double[2];

x0[0] = -2;

x0[1] = 4;

double x00 = -2;

double x01 = 4;

double jd = 0.0001;//精度

bool biaozhi1 = true;//设置是否循环的标志

bool biaozhi2 = true;

double y;//定义关于x的函数值

double[] g = new double[2];//定义函数在点x处的梯度值

double[] p = new double[2];

double ggm, ggo=0;

double jixiao;

double x1=0, x2=0, fy=0;

int i;

for (int j = 0; j < 2; j++)

{

x[j] = x0[j];

}

while (biaozhi1)

{

i = 0;

while (biaozhi2)

{

y = JS.f(x[0], x[1]);//调用目标函数

g[0] = JS.g0(x[0], x[1]);//调用梯度函数,求在点x处的梯度值 g[1] = JS.g1(x[0], x[1]);

ggm = g[0] * g[0] + g[1] * g[1];

if (ggm <= jd)

{

x1 = x[0];

x2 = x[1];

fy = y;

biaozhi1 = false;

break;

}

else

{

if (i == 0)

{

p[0] = -g[0];

p[1] = -g[1];

}

else

{

p[0] = -g[0] + (ggm / ggo) * p[0];

p[1] = -g[1] + (ggm / ggo) * p[1];

}

//调用0.618法子函数,找出极小点

jixiao = JS.JiXiao1(x[0], x[1], p[0], p[1], x00, x01); x[0] = x[0] + jixiao * p[0];

x[1] = x[1] + jixiao * p[1];

ggo = ggm;

i++;

if (i >= 2)

{

break;

}

else

{

continue;

}

}

}

}

txtx1.Text = x1.ToString();

txtx2.Text = x2.ToString();

txty.Text = fy.ToString();

}

}

②子函数程序:

using System;

using System.Data;

using System.Configuration;

using System.Web;

using System.Web.Security;

using System.Web.UI;

using System.Web.UI.WebControls;

using System.Web.UI.WebControls.WebParts;

using System.Web.UI.HtmlControls;

using System.Collections;

///

/// JiSuan 的摘要说明

///

public class JiSuan

{

public JiSuan()

{

//

// TODO: 在此处添加构造函数逻辑

//

}

//目标函数

public double f(double z1,double z2)

{

double zz = 1.5 * z1 * z1 + 0.5 * z2 * z2 - z1 * z2 - 2 * z1;

return zz;

}

//求梯度值函数

public double g0(double z1, double z2)

{

double zz = 3 * z1 - z2 - 2;

return zz;

}

public double g1(double z1, double z2)

{

double zz = z2 - z1;

return zz;

}

//倍增半减函数

public ArrayList SuoJian(double x0, double t0, double p, double q, double xx0, double xx1, double p0, double p1,double x00,double x01)

{

double t1;

double f0, f1;

double x1, t2;

double f2;

double a = 0, b = 0, c = 0;

double fa, fc, fb;

double t3, f3;

bool biaozhi1 = true;//设置是否循环的标志

bool biaozhi2 = true;

double x_0t1, x_1t1, x_0t2, x_1t2;

x_0t1 = xx0 + t1 * p0;

x_1t1 = xx1 + t1 * p1;

f0 = f(x00, x01);//调用目标函数

f1 = f(x_0t1, x_1t1);

if (f0 > f1)

{

while (biaozhi1)

{

x1 = p * x0;

t2 = t1 + x1;

x_0t2 = x_0t1 + t2 * p0;

x_1t2 = x_1t1 + t2 * p1;

f2 = f(x_0t2, x_1t2);//调用目标函数

if (f1 > f2)

{

t0 = t1;

t1 = t2;

f0 = f1;

f1 = f2;

continue;

}

else

{

a = t0;

c = t1;

b = t2;

fa = f0;

fc = f1;

fb = f2;

break;

}

}

}

else

{

t3 = t1;

f3 = f1;

t1 = t0;

f1 = f0;

t0 = t3;

while (biaozhi2)

{

x1 = q * x0;

t2 = t1 - x1;

x_0t2 = x_0t1 + t2 * p0;

x_1t2 = x_1t1 + t2 * p1;

f2 = f(x_0t2, x_1t2);//调用目标函数

if (f1 > f2)

{

t0 = t1;

t1 = t2;

f0 = f1;

f1 = f2;

continue;

}

else

{

a = t2;

c = t1;

b = t0;

fa = f2;

fc = f1;

fb = f0;

break;

}

}

}

ArrayList list = new ArrayList();

list.Add(a);

list.Add(b);

return list;

}

//0.618法函数

public double JiXiao1(double x0, double x1, double p0, double p1,double x00,double x01)

{

double jd = 0.0001;//精度

double p = 0.618;

double t1, t2;

double jixiao = 0;//极小点

bool biaozhi1 = true;//设置是否循环标志

bool biaozhi2 = true;

double xx0 = 0.01;//步长

double tt0 = 0, pp = 2, qq = 0.5;

double x_0t1, x_1t1, x_0t2, x_1t2;

ArrayList list1 = new ArrayList();

//调用倍增半减函数获取极小区间

list1 = (ArrayList)SuoJian(xx0, tt0, pp, qq, x0, x1, p0, p1, x00, x01);

double a = double.Parse(list1[0].ToString());

double b = double.Parse(list1[1].ToString());

while (biaozhi1)

{

t1 = a + (1 - p) * (b - a);

t2 = a + p * (b - a);

x_0t1 = x00 + t1 * p0;

x_1t1 = x01 + t1 * p1;

x_0t2 = x_0t1 + t2 * p0;

x_1t2 = x_1t1 + t2 * p1;

f1 = f(x_0t1, x_1t1);//调用目标函数

f2 = f(x_0t2, x_1t2);

while (biaozhi2)

{

if (f1 < f2)

{

b = t2;

t2 = t1;

f2 = f1;

t1 = a + (1 - p) * (b - a);

x_0t1 = x00 + t1 * p0;

x_1t1 = x01 + t1 * p1;

f1 = f(x_0t1, x_1t1);//调用目标函数

if (Math.Abs(b - a) > jd)

{

continue;

}

else

{

biaozhi1 = false;

break;

}

}

else if (f1 == f2)

{

a = t1;

b = t2;

if (Math.Abs(b - a) > jd)

{

break;

}

else

{

biaozhi1 = false;

break;

}

}

else if (f1 > f2)

{

a = t1;

t1 = t2;

f1 = f2;

t2 = a + p * (b - a);

x_0t2 = x_0t1 + t2 * p0;

x_1t2 = x_1t1 + t2 * p1;

f2 = f(x_0t2, x_1t2);//调用目标函数

if (Math.Abs(b - a) > jd)

{

continue;

}

else

{

biaozhi1 = false;

break;

}

}

}

}

jixiao = (a + b) / 2;

return jixiao;

}

}

三、程序运行界面及结果

(1)0.618法和抛物线法求ψ(t)=sint,初始点t0=1

(2)共轭梯度法求函数极小点:f(x)=1.5x12+0.5x22-x1x2-2x1;初始点为(-2,4)

北航最优化方法大作业参考

北航最优化方法大作业参考

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

数学实验“线性方程组的最速下降法与共轭梯度法解法”实验报告(内含matlab程序代码)

西京学院数学软件实验任务书

实验五实验报告 一、实验名称:最速下降法与共轭梯度法解线性方程组。 二、实验目的:进一步熟悉理解掌握最速下降法与共轭梯度法解法思路,提高matlab 编程能力。 三、实验要求:已知线性方程矩阵,应用最速下降与共轭梯度法在相关软件编程求解线性方程组的解。 四、实验原理: 1.最速下降法: 从某个初始点)0(X 出发,沿)(X f 在点)0(X 处的负梯度方向 )0()0()0()(AX b X f r -=-?= 求得)(X f 的极小值点)1(X , 即 )(min )0()0(0 r X f λλ+> 然后从)1(X 出发,重复上面的过程得到)2(X 。如此下去,得到序列{)(k X } )(...)()()()1()0(k X f X f X f >>> 可以证明,从任一初始点)0(X 出发, 用最速下降法所得到的序列{)(k X }均收敛于问题使X 最小化)(X f 的解,也就是方程组b AX =的解。其收敛速度取决于 1 1 λλλλ+-n n ,其中1λ ,n λ分别

为A 的最小,最大特征值。最速下降法迭代格式:给定初值)0(X , )(k X 按如下方法决定: ()) ()(1)(k )()()()(k ) ()(X ,,)(k k k k T k k T k k k k r X Ar r r r AX b X f r λλ+=> <><=-=-?=+ 2.共轭梯度法 其基本步骤是在点)(k X 处选取搜索方向)(k d , 使其与前一次的搜索方向)1(-k d 关于A 共轭,即 (1)()(1),0k k k d d Ad --<>= 然后从点)(k X 出发,沿方向)(k d 求得)(X f 的极小值点 )1(+k X , 即 )(min )() ()(0 )1(k d X f X f k k λλ+=>+ 如此下去, 得到序列{)(k X }。不难求得0,)1()(>=<-k k Ad d 的解为 ) () 1()1()()() () 1(,,k k k k k k k d Ad d d AX b X X > <>-<+=--+ 注意到)(k d 的选取不唯一,我们可取

最优化方法大作业答案

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-- 选取最优解:

最优化方法实验报告(2)

最优化方法实验报告Numerical Linear Algebra And Its Applications 学生所在学院:理学院 学生所在班级:计算数学10-1 学生姓名:甘纯 指导教师:单锐 教务处 2013年5月

实验三 实验名称: 无约束最优化方法的MATLAB 实现 实验时间: 2013年05月10日 星期三 实验成绩: 一、实验目的: 通过本次实验的学习,进一步熟悉掌握使用MATLAB 软件,并能利用该软件进行无约束最优化方法的计算。 二、实验背景: (一)最速下降法 1、算法原理 最速下降法的搜索方向是目标函数的负梯度方向,最速下降法从目标函数的负梯度方向一直前进,直到到达目标函数的最低点。 2、算法步骤 用最速下降法求无约束问题n R x x f ∈,)(min 的算法步骤如下: a )给定初始点)0(x ,精度0>ε,并令k=0; b )计算搜索方向)()()(k k x f v -?=,其中)()(k x f ?表示函数)(x f 在点)(k x 处的梯度; c )若ε≤)(k v ,则停止计算;否则,从)(k x 出发,沿)(k v 进行一维搜索, 即求k λ,使得)(min )()()(0 )()(k k k k v x f v x f λλλ+=+≥; d )令1,)()()1(+=+=+k k v x x k k k k λ,转b )。

(二)牛顿法 1、算法原理 牛顿法是基于多元函数的泰勒展开而来的,它将 )()]([-)(1)(2k k x f x f ??-作为搜索方向,因此它的迭代公式可直接写出 来: )()]([)(1)(2)()(k k k k x f x f x x ??-=- 2、算法步骤 用牛顿法求无约束问题n R x x f ∈),(min 的算法步骤如下: a )给定初始点)0(x ,精度0>ε,并令k=0; b )若ε≤?)()(k x f ,停止,极小点为)(k x ,否则转 c ); c )计算)()]([,)]([)(1)(2)(1)(2k k k k x f x f p x f ??-=?--令; d )令1,)()()1(+=+=+k k p x x k k k ,转b )。 (三)共轭梯度法 1、算法原理 共轭梯度法是利用目标函数梯度逐步产生共轭方向作为线搜索方向的方法,每次搜索方向都是在目标函数梯度的共轭方向,搜索步长通过一维极值算法确定。 2、算法步骤 a )给定初始点)0(x ,精度0>ε;

共轭梯度法C语言(西安交大)

#include #include #define N 10 /*定义矩阵阶数*/ void main() { int i,j,m,A[N][N],B[N]; double X[N],akv[N],dka[N],rk[N],dk[N],pk,pkk,ak,bk; for(i=0;i

printf("\n"); printf("input the Maxtr X\n"); /*给X0输入一个初始向量*/ for(i=0;i

最优化方法大作业

发动机空燃比控制器 引言:我主要从事自动化相关研究。这里介绍我曾经接触过的发动机空燃比控制器设计中的优化问题。 发动机空燃比控制器设计中的最优化问题 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 的差异尽可能的小。所以原问

最优化原理大作业

基于粒子群算法的神经网络在电液伺服系统中的应用 摘要:由于人工神经网络在解决具有非线性、不确定性等系统的控制问题上具有极大的潜力,因而在控制领域正引起人们的极大关注,并且已在一些响应较慢的过程控制中获得成功应用。由于电液伺服系统属 于非线性系统,因此本文利用神经网络控制电液伺服系统,并利用粒子群优化算法训练该神经网络的 权值。通过对神经网络的优化实现对电液伺服系统的控制。 关键词:神经网络电液伺服系统粒子群算法优化 近年来,由于神经网络具有大规模并行性、冗余性、容错性、本质的非线性及自组织自学习自适应能力,所以已成功地应用于众多领域。但在具有复杂非线性特性的机电设备的实时控制方面,虽然也有一些神经网络技术的应用研究,但距实用仍有一段距离。电液伺服系统就属于这类设备[1]。 神经网路在用于实时控制时,主要是利用了网络所具有的其输人——输出间的非线性映射能力。它实际上是通过学习来逼近控制对象的动、静态特性。也就是构造实际系统的神经网络模型[2]。本文利用神经网络控制一电液伺服系统,并利用粒子群优化算法训练该神经网络的权值,将结果与BP神经网络控制该系统的结果进行比较。从而得在电液伺服系统中引入神经网络是可行的。 1、粒子群算法 粒子群优化算法(Particle Swarm optimization, PSO)是一种进化计算技术, 由Eberhart博士和kennedy博士发明, 源于对鸟群捕食的行为研究, 粒子群优化算法的基本思想是通过群体中个体之间的协作和信息共享来寻找最优解[3]。算法最初受到飞鸟和鱼类集群活动的规律性启发,利用群体智能建立了一个简化模型,用组织社会行为代替了进化算法的自然选择机制,通过种群间个体协作来实现对问题最优解的搜索[4]。 在找到这两个最优值时, 粒子根据如下的公式来更新自己的速度和新的位置 v[]=v[]+c1*rand()*(pbest[]-present[]) + c2*rand()*(gbest[]-present[]) present[]=persent[]+v[] 式中ω为惯性权重,ω取大值可使算法具有较强的全局搜索能力,ω取小值则算法倾向于局部搜索。一般的做法是将ω初始取0.9并使其随迭代次数的增加而线性递减至0.4,这样就可以先侧重于全局搜索,使搜索空间快速收敛于某一区域,然后采用局部精细搜索以获得高精度的解;c1、c2为两个学习因子,一般取为2;randl和rand2为两个均匀分布在(0,l)之间的随机数;i=1,2,?,m;k=1,2,?,d。另外,粒子在每一维的速度Vi都被一个最大速度Vmax所限制。如果当前粒子的加速度导致它在某一维的速度 超过该维上的最大速度Vmax,则该维的速度被限制为最大速度[5]。 粒子群算法流程如下: (一)初始化粒子群。设群体规模为m,在允许的范围内随机设置粒子的初始位置和速 度。 (二)评价每个粒子的适应值。 (三)调整每一个粒子的位置和速度。 (四)如果达到最大迭代次数genmax或误差达到最初设定数值终止迭代,否则返回(2)。 2、神经网络 神经网络一般由输入层、隐含层、输出层组成。对于输入信号,先向前传播到隐节点,经过节点作用函数后,再把隐节点的输出信息传播到输出节点,最后输出结果。节点的作用函数通常选取S 型函数f(x)=1/(1+e-x)。神经网络算法的学习过程分为正

最优化课程设计--共轭梯度法算法分析与实现

最优化课程设计--共轭梯度法算法分析与实现(设计程序) 题目共轭梯度法算法分析与实现 班级 / 学号 14140101/2011041401011 学生姓名黄中武指导教师王吉波王微微 课程设计任务书 课程名称最优化方法课程设计院(系) 理学院专业信息与计算科学 课程设计题目共轭梯度法算法分析与实现课程设计时间: 2014 年 6月 16日至 2014 年 6月 27日 课程设计的要求及内容: [要求] 1. 学习态度要认真,要积极参与课程设计,锻炼独立思考能力; 2. 严格遵守上机时间安排; 3. 按照MATLAB编程训练的任务要求来编写程序; 4. 根据任务书来完成课程设计论文; 5. 报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”; 6. 报告上交时间:课程设计结束时上交报告; 7. 严禁抄袭行为,一旦发现,课程设计成绩为不及格。 一、运用共轭梯度法求解无约束最优化问题 要求:1)了解求解无约束最优化问题的共轭梯度法; 2)绘出程序流程图; 3)编写求解无约束最优化问题的共轭梯度法MATLAB程序; 4)利用编写文件求解某无约束最优化问题;

5)给出程序注释。 指导教师年月日 负责教师年月日 学生签字年月日 沈阳航空航天大学 课程设计成绩评定单 课程名称最优化理论与算法课程设计院(系) 理学院专业信息与计算科学课程设计题目共轭梯度法算法分析与实现学号 2011041401011 姓名黄中武指导教师评语: 课程设计成绩 指导教师签字 年月日 最优化方法课程设计沈阳航空航天大学课程设计用纸目录 目录 一、正 文 (1) 二、总结 ............................................................... 8 参考文 献 ............................................................... 9 附录 .. (10) 第 I 页 最优化方法课程设计沈阳航空航天大学课程设计用纸正文 一、正文 一无约束最优化问题的共轭梯度法

北航惯性导航大作业

惯性导航基础课程大作业报告(一)光纤陀螺误差建模与分析 班级: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;

实验2 最速下降法和共轭梯度法的程序设计

实验2 最速下降法和共轭梯度法的程序设计 一、实验目的 1、熟悉无约束优化问题的最速下降算法和共轭梯度法。 2、培养matlab 编程与上机调试能力。 二、实验课时:2个课时 三、实验准备 1、预习无约束优化问题的最速下降算法和共轭梯度法。 2、熟悉matlab 软件的基本操作及程序编写。 四、实验内容 课堂实验演示 根据最速下降法编写程序,求函数 21222121342)(min x x x x x x x f -++-= 的极小值,其中初始点为()01,1T x = 算法步骤如下: Step1::给出初始点0x ,和精度1;0k ε<<=; Step2:计算()k f x ?,如果()k f x ε?≤,则停止迭代,输出结果;否则转step3; Step3:令下降方向()k k d f x =-?,计算步长因子k λ使得0()min ()k k k k k f x d f x d λλλ≥+=+,令1,1k k k k x x d k k λ+=+=+,转step2。 其程序如下: function [x,iter,val,dval] = Steepest_Descent_Method(x,eps) k = 1; dy = grad_obj(x); x_mat(:,1) = x;%存储每一次迭代得到的点x while norm(dy)>eps d = -dy; % 搜索方向 lambda = line_search(x,d);%步长 x = x + d*lambda; k = k + 1; x_mat(:,k) = x; dy = grad_obj(x); end iter = k - 1; val = obj(x);%目标函数在极值点处的函数值

作业4-FR共轭梯度法

最优化方法第四次作业 题目:利用FR-共轭梯度法求解无约束优化问题222 12122min ()44412x R f x x x x x x ∈=+--。初始点(0)(0.5,1).T x =- () ()T k k T k k k k k k k g g g g k d g k g d 111110.0,;0,-----=???≥+-=-=ββ 一、程序 function [x,val,k]=frcg(fun,gfun,x0) %功能:用FR 共轭梯度法求解无约束问题min f (x ) %输入:x0是初始点,fun,gfun 分别是求目标函数和梯度 %输出:x,val 分别是近似最优点和最优值,k 是迭代次数 maxk=5000; rho=0.6; sigma=0.4; k=0; epsilon=1e-4; n=length(x0); while (k=0.0) d=-g; end end if (norm(d)

while (m<20) %用Armijo 搜索求步长 if (feval(fun,x0+rho^m*d)> x0=[-0.5,1]'; >> [x,val,k]=frcg('fun','gfun',x0) x = 1.0000 2.0000 val = -12.0000 k = 10 即22212122min ()44412x R f x x x x x x ∈=+--的极小值点x=[1;2];minf(x)= -12。

最优化方法大作业答案

武工院你们懂的 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

列成表格: 00001216 100114 60105122001112----- 可见此表已具备1°,2°,3°三个特点,可采用单纯形法。首先从底行中选元素-1,由2/2,6/2,6/4最小者决定选第一行第一列的元素2,标以记号,迭代一次得 0000 1 2 121023 10 40116201002 1 21 211-------- 再从底行中选元素-2/3,和第二列正元素1/2,迭代一次得 1 002 1232 30210231 040116201002121211-- ----- 再从底行中选元素-3,和第二列正元素2,迭代一次得 4002 3 03410120280114042001112--- 再迭代一次得

10 23021 062 21023 1010 213 000421 2 10 13- - 选取最优解: 01=x 42=x 23=x 3. 试用DFP 变尺度法求解下列无约束优化问题。 min f (X )=4(x 1-5)2+(x 2-6)2 取初始点X=(8,9)T ,梯度精度ε=0.01。 解:取I H =0,初始点()T X 9,8= 2221)6()5(4)(-+-=x x x f ??????--=?122408)(21x x x f ???? ??=?624)() 0(x f T x f d )6,24()()0()0(--=-?= )0(0)0()1(d x x α+= T )69,248(00αα--= ])669()5248(4min[)(min 2020)0(0)0(--+--?=+αααd x f )6()63(2)24()2458(8) (00)0(0)0(=-?-+-?--=+ααααd d x df 13077.013017 0≈= α ???? ??=???? ??--?+???? ??=21538.886153.462413077.098)1(x

北航数值分析大作业第二题精解

目标:使用带双步位移的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 #include #include #define N 10 #define E 1.0e-12 #define MAX 10000 //以下是符号函数 double sgn(double a) { double z; if(a>E) z=1; else z=-1; return z; } //以下是矩阵的拟三角分解 void nishangsanjiaodiv(double A[N][N]) { int i,j,k; int m=0; double d,c,h,t; double u[N],p[N],q[N],w[N]; for(i=0;i

数值分析实验报告1——Hilbert矩阵的求解

数值分析课程实验报告 题目:病态线性方程组的求解 理论分析表明,数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解 Hx = b ,期中H 是Hilbert 矩阵,()ij n n H h ?=,1 1 ij h i j =+-,i ,j = 1,2,…,n 1. 估计矩阵的2条件数和阶数的关系 2. 对不同的n ,取(1,1,,1)n x =∈ ,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,比较结果。 3. 结合计算结果,试讨论病态线性方程组的求解。 解答过程 1.估计矩阵的2-条件数和阶数的关系 矩阵的2-条件数定义为:1 222 ()Cond A A A -=?,将Hilbert 矩阵带入有: 1222 ()Cond H H H -=? 调用自编的Hilbert_Cond 函数对其进行计算,取阶数n= 50,可得从1阶到50阶的2-条件数,以五位有效数字输出,其中前10项见表1。 表1.前十阶Hilbert 矩阵的2-条件数 从表1可以看出,随着阶数每递增1,Hilbert 矩阵的2-条件数都至少增加一个数量级,但难以观察出明显的相依规律。故考虑将这些数据点绘制在以n 为横轴、Cond (H )2为纵轴的对数坐标系中(编程用Hilbert_Cond 函数同时完成了这个功能),生成结果如图1。

图1.不同阶数下Hilbert矩阵的2-条件数分布 由图可见,当维数较小时,在y-对数坐标系中Cond(H)2与n有良好的线性关系;但n超过10后,线性趋势开始波动,n超过14后更是几乎一直趋于平稳。事实上,从n = 12开始,系统便已经开始提出警告:“Warning: Matrix is close to singular or badly scaled.Results may be inaccurate.”。也就是说,当n较大时,H矩阵已经接近奇异,计算结果可能是不准确的。通过查阅相关资料,我找到了造成这种现象的原因:在matlab中,用inv函数求条件数过大的矩阵的逆矩阵将是不可靠的。而调用系统自带的专门对Hilbert矩阵求逆的invhilb(n)函数则不存在这个问题,生成结果如图2。 图2. 修正后的不同阶数下Hilbert矩阵的2-条件数分布

肌组织实验报告

竭诚为您提供优质文档/双击可除 肌组织实验报告 篇一:表面肌实验报告 武汉理工大学 现代数字信号处理在前沿学科中的应用实验报告 基于semg时域特征的动作识别 学院:信息工程学院 学号:姓名: 班级:电子154 实验基于semg时域特征特的动作识别 一、实验目的 1.了解肌电信号常用的时域分析方法; 2.利用mATLAb对肌电信号进行去噪、特征提取及动作识别; 二、实验设备 1.wi-Fi表面肌电信号采集卡; 2.32位windowsxp台式机(matlab7.0软件); 3.802.11b/g无线网卡;

三、实验内容 (1)学习信号的基本去噪方法,并用mATLAb实现; (2)学习肌电信号常用的时域特征并利用matlab来进行波形长度(wL)符号改变数(ssc)、过零点(Zc)、威尔 逊赋值(wAmp)等特征的提取; (3)学习神经网络信号处理方法,掌握bp神经网络的用法,将其用于肌电信号的动作识别。 学习以上三个部分,最终完成一整套肌电信号去噪、特征提取(选取一种特征)、基于特征的动作识别的mATLAb程序。 四、实验原理 (1)小波去噪 小波去噪方法是一种建立在小波变换基础上的新兴算法,基本思想是根据噪声在不同频带上的小波分解系数具有不同强度分布的特点,将各频带上的噪声对应的小系数去除,保留原始信号的小波分解系数,然后对处理后系数进行小波重构,得到纯净信号。 小波去噪的基本原理图如下 (2)特征提取 时域分析是将肌电信号看成均值为零,而方差随着信号强度的变化而变化的随机信号。时域特征的计算复杂度低,提取比较方便。

最常用的方法有:方差,过零点数(Zerocrossing,Zc),willison幅值(willisonAmplitude,wAmp),绝对值平均值(meanAbsoluteValue,mAV)和波形长度(wavelength,wL)等。在实际应用中,为了让特征可以包含更多的信息,往往选择用不同的时域特征组合形成联合特征向量。我们主要介绍一下几种方法: 过零率(Zc):为波形通过零线的次数,从一定程度上反映了信号的频率特性。为了降低零点引入的噪声,往往会引入一个阈值δ。计算方式如下: sgn(?xk?xk?1),(xk?xk?1??)(1)willison幅值:是由willison提出一种对表面肌电信号的幅值变化数量进行计 算的方法,经过后人的研究,对willison幅值的阈值有了明确的范围限定,目前认为50~100?V是最合适的阈值范围。其数学表示公式如公式(3-3)。 wAmp??fxi?xi?1 t?1n(2) ?1f(x)???0其中:ifx?阈值otherwise 波形长度(wL):它是对某一分析窗中的波形长度的统计,波长可以体现该样本的持续时间、幅值、频率的特征。 1n?1 wL??x(i?1)?x(i)ni?1(3)符号改变斜率(ssc):为信号的的频率性能提供了一些附加信息,对于3个连续的采样

大连理工优化方法大作业MATLAB编程

function [x,dk,k]=fjqx(x,s) flag=0; a=0; b=0; k=0; d=1; while(flag==0) [p,q]=getpq(x,d,s); if (p<0) b=d; d=(d+a)/2; end if(p>=0)&&(q>=0) dk=d; x=x+d*s; flag=1; end k=k+1;

if(p>=0)&&(q<0) a=d; d=min{2*d,(d+b)/2}; end end %定义求函数值的函数fun,当输入为x0=(x1,x2)时,输出为f function f=fun(x) f=(x(2)-x(1)^2)^2+(1-x(1))^2; function gf=gfun(x) gf=[-4*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),2*(x(2)-x(1)^2)]; function [p,q]=getpq(x,d,s) p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s'; q=gfun(x+d*s)*s'-0.60*gfun(x)*s'; 结果: x=[0,1]; s=[-1,1]; [x,dk,k]=fjqx(x,s) x =-0.0000 1.0000 dk =1.1102e-016 k =54

function f= fun( X ) %所求问题目标函数 f=X(1)^2-2*X(1)*X(2)+2*X(2)^2+X(3)^2+ X(4)^2- X(2)*X(3)+2*X(1)+3*X(2)-X(3); end function g= gfun( X ) %所求问题目标函数梯度 g=[2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X(3)-X(2)-1,2*X(4)]; end function [ x,val,k ] = frcg( fun,gfun,x0 ) %功能:用FR共轭梯度法求无约束问题最小值 %输入:x0是初始点,fun和gfun分别是目标函数和梯度 %输出:x、val分别是最优点和最优值,k是迭代次数 maxk=5000;%最大迭代次数 rho=0.5;sigma=0.4;

共轭梯度法程序

一、共轭梯度法 共轭梯度法(Conjugate Gradient)是共轭方向法的一种,因为在该方向法中每一个共轭向量都是依靠赖于迭代点处的负梯度而构造出来的,所以称为共轭梯度法。由于此法最先由Fletcher和Reeves (1964)提出了解非线性最优化问题的,因而又称为FR 共轭梯度法。由于共轭梯度法不需要矩阵存储,且有较快的收敛速度和二次终止性等优点,现在共轭梯度法已经广泛地应用于实际问题中。共轭梯度法是一个典型的共轭方向法,它的每一个搜索方向是互相共轭的,而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合,因此,存储量少,计算方便,效果好。 二、共轭梯度法的原理 设有目标函数 f(X)=1/2X T HX+b T X+c 式1 式中,H作为f(X)的二阶导数矩阵,b为常数矢量,b=[b1,b2,b3,...b n]T 在第k次迭代计算中,从点X(k)出发,沿负梯度方向作一维搜索,得 S(K)=-?f(X(k))式2 X(k+1)=X(k)+ɑ(k)S(k) 式3 在式中,ɑ(k)为最优步长。 设与S(k)共轭的下一个方向S(k+1)由点S(k)和点X(k+1)负梯度的线性组

合构,即 S (k+1)=-?f (X (k+1))+β(k)S (k) 式4 根据共轭的条件有 [S (k)]T ?2f (X (k))S (k+1)=0 式5 把式2和式4带入式5,得 -[?f(X (k))]T ?2f (X (k))[-?f (X (k+1))+β(k)S (k) ]=0 式6 对于式1,则在点X (k)和点X (k+1)的梯度可写为 ?f(X (k))=HX (k)+b 式7 ?f (X (k+1))=HX (k+1)+b 式8 把上面两式相减并将式3代入得 ɑ(k)H S (k)=?f (X (k+1))-?f(X (k)) 式9 将式4和式9两边分别相乘,并代入式5得 -[?f (X (k+1))+β(k)?f(X (k))]T [?f (X (k+1))-?f(X (k)]=0 式10 将式10展开,并注意到相邻两点梯度间的正交关系,整理后得 β (k ) =2 2 ||))((||||))1((||k X f k X f ?+? 式11 把式11代入式4和式3,得 S (k+1)=-?f (X (k))+β (k ) S (k ) X (k+1)=X (k )+ɑ(k )S (k ) 由上可见,只要利用相邻两点的梯度就可以构造一个共轭方向。以这种方式产生共轭方向并进行迭代运算的方法,即共轭梯度法。

结构优化设计大作业(北航)

《结构优化设计》 大作业报告 实验名称: 拓扑优化计算与分析 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插值模型。

相关文档
最新文档