机械优化设计二次插值法流程图
机械优化设计-第三章一维优化方法
机械优化设计
• 第四次缩小区间： 第四次缩小区间： • 令 x2=x1=0.764, , f2=f1=0.282 • x1=0.472+0.382*(0.944-0.472)=0.652, f1=0.223 • 由于f1<f2, 故新区间 由于f 故新区间[a,b]=[a, x2]=[0.472, 0.764] • 因为 b-a=0.764-0.472=0.292>0.2, 应继续缩小区间。 , 应继续缩小区间。 第五次缩小区间： 第五次缩小区间： f2=f1=0.223 令 x2=x1=0.652, x1=0.472+0.382*(0.764-0.472)=0.584, f1=0.262 由于f 故新区间[a,b]=[x1,b]=[0.584, 0.764] 由于f1>f2, 故新区间 因为 b-a=0.764-0.584=0.18<0.2, 停止迭代。 程序演示 , 停止迭代。 极小点与极小值： 极小点与极小值： x*=0.5*(0.584+0.764)=0.674,
x2 = a + 0.618(b − a), y2 = f ( x2 )
f
b = x2 , x2 = x1, y2 = y1
x1 = a + 0.382(b − a), y1 = f ( x1 )
y1 < y2
否
是
y1 y2
x
b
a = x1 , x1 = x2 , y1 = y2
x2 = a + 0.618(b − a), y2 = f ( x2 )
7
机械优化设计
h0
x2
机械优化设计
2.前进搜索 加大步长 h＝2 h ，产生新点x3= x2+ 2h0 ； （a）如y2<y3，则函数在[x1，x3]内 必有极小点，令a= x1，b= x3搜索 区间为[a，b] ； （b）如y2>y3， 令x1=x2 ，y1=y2 ； x2=x3 ，y2=y3 ； h=2h 重新构造新点x3=x2+h，并比较y2、 y3的大小，直到y2<y3。
机械优化设计(全,120张)
| F ( X (k1) ) F ( X (k) ) |
或
|
F(
X
(k1) ) F ( X | F(X (k)) |
(k)
)
|
(三)梯度准则
目标函数在迭代点的梯度模已达到充分小， 即：
|| F ( X (k1) ||
式中的ε为结定的迭代精度或允许误差，可 根据设计要求预先给定。
1.3 优化设计问题的基本解法
求解优化问题的方法：
解析法
数学模型复杂时不便求解
数值法
可以处理复杂函数及没有数学表达式 的优化设计问题
寻求极值点的搜索过程
1.4、迭代计算的终止准则
在迭代过程中，如果根据一个迭代公式能够计 算出接近精确解的近似解，即近似解序列X（k) 有极限
lim X (k) X * k
在设计这根主轴时， 有两个重要因素需要 考虑。一是主轴的自 重；一是主轴伸出端 c点的挠度。
1.11、设计变量
在优化设计的过程中，不断进行修改、调整， 一直处于变化的参数称为设计变量。
设计变量的全体实际上是一组变量，可用一 个列向量表示：
x x1 x2 ... xn T
本问题中，当主轴的材料选定时，其设计方案 由四个设计变量决定。孔径d、外径D、跨距l及 外伸端长度a。由于机床主轴内孔用于通过待加 工的棒料，其大小由机床型号决定。不作为设 计变量。故设计变量取为：
2、现代设计方法
现代设计方法是以电子计算机 为手段,运用工程设计的新理论和 新方法,使计算结果达到最优化,使 设计过程实现高效化和自动化.
运用现代设计方法可以适应市 场剧烈竞争的需要,提高设计质量 和缩短设计周期.
二、优化设计方法
机械优化设计 多项式插值和分段插值1
一、插值思想简介
5/27
一、插值思想简介
引例1、函数查表问题
求标准正态分布函数值F(2.3456789) 由标准正态分布函数值表可以得到： F(2.34)=0.99036; F(2.35)=0.99061 2.3456789接近2.35 故F(2.3456789)约等于0.99061
在对精度要求较高时，这种处理方法受到质疑 问题：利用一个表格给出的函数值，计算未给出的函数值
设函数f(x)在n+1个节点x0,x1,…,xn处的函数值已知， 为 y0,y1,…,yn 。 要求一个分段( 共 n段)线性函数q(x)，使其满足：
q(xi)=yi , i=0,1,…,n.
17/27
二、一维插值方法
这n+1个点 x0,x1,…,xn
称为节点
y
o x0
xj-1 xj xj+1 xn
n=3时
axis([-6,6,-1.5 2]); n=9时
高阶振荡
n=15时
15/27
二、一维插值方法
2、分段多项式插值 由于高次插值多项式的振荡缺陷，促使人们转而寻求简单的低 次多项式插值。分段插值的基本思想是将插值区间划分为若干 个小区间，然后在每个小区间上做满足一定条件的低阶插值。 2.1分段线性插值
Li
(x)
(x x0 )( x x1)(x (xi x0 )( xi x1)(xi
xi1)( x xi1)(x xn ) xi1)( xi xi1)( xi xn
)
易知n次多项
式Li
(
x)满足Li
(
x
j
)
1，i 0，i
线性插值与二次插值公式ppt课件
MATLAB计算程序
1
x=0:.6:1.8; y=erf(x);
0 .8
x=x';y=y';
A=[ones(4,1) x x.^2 x.^3]; 0.6
p=A\y;
0 .4
a0=p(1);a1=p(2); 0 .2
a2=p(3);a3=p(4);
t=0:.2:2;
0
0
0 .5
1
1 .5
2
u=a0+a1*t+a2*t.^2+a3*t.^3;
plot(x,y,'o',t,u)
12
拉格朗日插值的基函数构造法
n=1 线性插值问题 x
x0
x1
已知函数表 f(x)
y0
y1
求满足: L1(x0)=y0 , L1(x1)=y1的线性插值多项式 L1(x)
由过两点直线方程,得
L1( x)
y0
y1 x1
y0 x0
(x
x0 )
化为等价形式
L1( x)
当 x∈(0.5, 1)时
Erf ( x) 1 [( x 0.5) 0.8427 (1 x) 0.5205] 1 0.5
当 x∈(1, 1.5)时
Erf ( x) 1 [( x 1) 0.9661 (1.5 x) 0.8427] 1.5 1
3
实际问题中遇到的函数f(x)有的表达式复杂，有 的只提供了离散点上的函数值或导数值。为了进 一步分析问题的性质和变化规律，自然希望找到 一种简单函数p(x)，能近似描述函数f(x)的变化规 律，又便于处理。把这个函数p(x)称作f(x)的近似 函数。
机械优化设计黄金分割法 外推法
郑州大学机械优化设计部分程序1.外推法2.黄金分割法3.二次插值法4.坐标轮换法5.随机方向法6.四杆机构优化设计1.外推法源程序：#include<>#include<>#define Rdouble fun(double x){ double m;m=x*x-10*x+36;return m;}void main(){double h0=R,y1,y2,y3,x1,x2,x3,h; x1=0;h=h0;x2=h;y1=fun(x1);y2=fun(x2);if(y2>y1){h=-h;x3=x1;y3=y1;x1=x2;y1=y2;x2=x3;y2=y3;}x3=x2+h;y3=fun(x3);while(y3<y2){h*=;x1=x2;y1=y2;x2=x3;y2=y3;x3=x2+h;y3=fun(x3);}printf("fun(%f)=%f,fun(%f)=%f,fun( %f)=%f\n",x1,y1,x2,y2,x3,y3);}运行过程及结果：fun=,fun=,fun=2.黄金分割法源程序：#include<>#include<>#define f(x) x*x*x*x-5*x*x*x+4*x*x-6*x+60double hj(double *a,double *b,double e,int *n){double x1,x2,s;if(fabs((*b-*a)/(*b))<=e)s=f((*b+*a)/2);else{x1=**(*b-*a);x2=*a+*(*b-*a);if(f(x1)>f(x2))*a=x1;else*b=x2;*n=*n+1;s=hj(a,b,e,n);}return s;}void main() {double s,a,b,e,m;int n=0;printf("输入a，b值和精度e值\n"); scanf("%lf %lf %lf",&a,&b,&e);s=hj(&a,&b,e,&n);m=(a+b)/2;printf("a=%lf,b=%lf,s=%lf,m=%lf,n=%d\n ",a,b,s,m,n);}运行过程及结果：输入a，b值和精度e值-35a=,b=,s=,m=,n=213.二次插值法源程序：#include<>#include<>int main(void){doublea1,a2,a3,ap,y1,y2,y3,yp,c1,c2,m; double j[3];int i,h=1;void finding(double a[3]);finding(j);a1=j[0];a2=j[1];a3=j[2];m=;double f(double x);y1=f(a1);y2=f(a2);y3=f(a3);for(i=1;1>=1;i++){c1=(y3-y1)/(a3-a1);c2=((y2-y1)/(a2-a1)-c1)/(a2-a3); ap=*(a1+a3-c1/c2);yp=f(ap);if(fabs((y2-yp)/y2)<m)break;else if((ap-a2)*h>0){ if(y2>=yp){a1=a2;y1=y2;a2=ap;y2=yp;}else{a3=ap;y3=yp;}}else if(y2>=yp){a3=a2;y3=y2;a2=ap;y2=yp;}else{a1=ap;y1=yp;} }double x,y;if(y2<=yp){x=a2;y=y2;}else{x=ap;y=yp;}printf("a*=%f\n",x); printf("y*=%f\n",y); return 0;}double f(double x) {double y;y=x*x-10*x+36;return y;}void finding(double a[3]){int h,i;double y[3];a[0]=0;h=1;a[1]=h;y[0]=f(a[0]);y[1]=f(a[1]); if(y[1]>y[0]){h=-h;a[2]=a[0];y[2]=y[0];do{a[0]=a[1];a[1]=a[2];y[0]=y[1];y[1]=y[2];a[2]=a[1]+h;y[2]=f(a[2]); h=2*h;}while(y[2]<y[1]);}else{for(i=1;i>=1;i++){a[2]=a[1]+h;y[2]=f(a[2]); if(y[2]>=y[1])break;h=2*h;a[0]=a[1];y[0]=y[1];a[1]=a[2];y[1]=y[2];}}return;}运行过程及结果：a*=y*=4.坐标轮换法源程序：#include <>#include <>#include <>float fun1(float x,float a,float b) {float y;y=x+a*b;return y;}float fun2(float x,float y){float z;z=4*(x-5)*(x-5)+(y-6)*(y-6);return z;}main(){floatd[100][3],x[100][3],xx[3],ax[100][3]; floata1,a2,a3,h,t,y1,y2,y3,e,a,b,l,fi;int i,k;printf("输入初始点坐标\n");scanf("%f%f",&x[0][1],&x[0][2]); e=;l=;x[2][1]=x[0][1];x[2][2]=x[0][2];k=0;k--;do{x[0][1]=x[2][1];x[0][2]=x[2][2];k++;for(i=1;i<=2;i++){if(i==1){d[i][1]=1; d[i][2]=0;}else{d[i][1]=0;d[i][2]=1;}h=;a1=0;a2=h;x[i][1]=fun1(x[i-1][1],d[i][1],a1); x[i][2]=fun1(x[i-1][2],d[i][2],a1); y1=fun2(x[i][1],x[i][2]);x[i][1]=fun1(x[i-1][1],d[i][1],a2); x[i][2]=fun1(x[i-1][2],d[i][2],a2); y2=fun2(x[i][1],x[i][2]);if(y2>y1){h=-h;a3=a1;y3=y1;a1=a2;a2=a3;y1=y2;y2=y3;}a3=a2+h;x[i][1]=fun1(x[i-1][1],d[i][1],a3);x[i][2]=fun1(x[i-1][2],d[i][2],a3); y3=fun2(x[i][1],x[i][2]);do{a1=a2;y1=y2;a2=a3;y2=y3;a3=a2+h;x[i][1]=fun1(x[i-1][1],d[i][1],a3);x[i][2]=fun1(x[i-1][2],d[i][2],a3); y3=fun2(x[i][1],x[i][2]);}while(y3<y2);for(;a1>a3;){t=a3;a3=a1;a1=t; t=y1;y3=y1;y1=t;}a=a1;b=a3;a1=b-l*(b-a);a2=a+l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a1); x[i][2]=fun1(x[i-1][2],d[i][2],a1); y1=fun2(x[i][1],x[i][2]);x[i][1]=fun1(x[i-1][1],d[i][1],a2); x[i][2]=fun1(x[i-1][2],d[i][2],a2); y2=fun2(x[i][1],x[i][2]);if(b<1e-3){for(;fabs(b-a)>e;){if(y1>=y2){a=a1;a1=a2;y1=y2;a2=a+l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a2);x[i][2]=fun1(x[i-1][2],d[i][2],a2);y2=fun2(x[i][1],x[i][2]);}else{b=a2;a2=a1;y2=y1;a1=b-l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a1);x[i][2]=fun1(x[i-1][2],d[i][2],a1);y1=fun2(x[i][1],x[i][2]);}}}else{for(;fabs((b-a)/b)>=e||fabs((y2-y1 )/y2)>=e;){if(y1>=y2){a=a1; a1=a2;y1=y2;a2=a+l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a2);x[i][2]=fun1(x[i-1][2],d[i][2],a2);y2=fun2(x[i][1],x[i][2]);}else{b=a2;a2=a1;y2=y1;a1=b-l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a1);x[i][2]=fun1(x[i-1][2],d[i][2],a1);y1=fun2(x[i][1],x[i][2]);}}}ax[k][i]=*(a+b);x[i][1]=fun1(x[i-1][1],d[i][1],ax[k][i]);x[i][2]=fun1(x[i-1][2],d[i][2],ax[k][i ]);}}while(sqrt(pow((x[2][1]-x[0][1]),2)+po w((x[2][2]-x[0][2]),2))>=1e-6);xx[1]=x[2][1];xx[2]=x[2][2];fi=fun2(xx[1],xx[2]);printf("最优解为\nx1*=%f\nx2*=%f\nf*=%f\nk=%d\n",xx[1],xx[2],fi,k); }运行过程及结果：输入初始点坐标89最优解为x1*=x2*=f*=k=25.随机方向法源程序：#include<>#include<>#include<>float f(float x,float y){float z;z=(x-2)*(x-2)+(y-1)*(y-1); return z;}float g1(float x,float y){float z;z=x*x-y;return z;}float g2(float x,float y){float z;z=x+y-2;return z;}void main(){int i,j;float k=8,c=,a0=-3,b0=3,a1=-3,b1=3; floatx[10],x0[10],xl[10],e[10],r[10],d[10], h,fl,f0,fx;while(g1(x0[0],x0[1])>0||g2(x0[0],x0[1 ])>0){x0[0]=a0+(rand()/*(b0-a0);x0[1]=a1+(rand()/*(b1-a1);}fl=f(x0[0],x0[1]);f0=f(x0[0],x0[1]);while(1){h=;j=1;r[0]=-1+(rand()/*(1-(-1));r[1]=-1+(rand()/*(1-(-1));e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);x[0]=x0[0]+h*e[0];x[1]=x0[1]+h*e[1];if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0){fx=f(x[0],x[1]);if(fx<fl){fl=fx;for(i=0;i<2;i++){d[i]=e[i];xl[i]=x[i];}}}while(j<=k){j++;r[0]=-1+(rand()/*(1-(-1));r[1]=-1+(rand()/*(1-(-1));e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);x[0]=x0[0]+h*e[0];x[1]=x0[1]+h*e[1];if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0) {fx=f(x[0],x[1]);if(fx<fl){fl=fx;for(i=0;i<2;i++){d[i]=e[i];xl[i]=x[i];}}}}x[0]=xl[0];x[1]=xl[1];while(1){h=*h;x[0]=x[0]+h*d[0];x[1]=x[1]+h*d[1];if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)break;fx=f(x[0],x[1]);if(fx<fl) fl=fx; else break;}do{x[0]=x[0]-h*d[0];x[1]=x[1]-h*d[1];h=*h;if(h<c)break;x[0]=x[0]+h*d[0];x[1]=x[1]+h*d[1];if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)continue;fx=f(x[0],x[1]);}while(fx>=fl);if(fabs((f0-fx)/f0)>=c){x0[0]=x[0];x0[1]=x[1];fl=fx;f0=fx;}elsebreak;}printf("输出最优解为\nx1*=%f,x2*=%f, y*=%f\n",x[0],x[1],fx);}运行过程及结果：输出最优解为x1*=,x2*=,y*=6.四杆机构优化设计源程序：#include<>#include<>#include<>#define Paiint g(float l1,float l2){if((-l1<=0)&&(-l2<=0)&&(6-l1-l2<=0)&&(1-l2-4<=0)&&(l2-l1-4<=0)&&(l1*l1+l2**l1*l2-16<=0)&&(36-l1*l1-l2**l1*l2<=0))return (1);elsereturn (0);}float fun(float x0[2]){floatf,a[31],b[31],r[31],p[31],q[31],w[31],x1 [2];int i;p[0]=acos(((1+x0[0])*(1+x0[0])-x0[ 1]*x0[1]+25)/(10+10*x0[0]));q[0]=acos(((1+x0[0])*(1+x0[0])-x0[ 1]*x0[1]-25)/(10*x0[1]));f=0;for(i=1;i<=30;i++){p[i]=p[0]+(Pai/60)*i;r[i]=sqrt(26-10*cos(p[i]));a[i]=acos((r[i]*r[i]+x0[1]*x0[1]-x 0[0]*x0[0])/(2*r[i]*x0[1]));b[i]=acos((r[i]*r[i]+24)/(10*r[i]) );q[i]=Pai-a[i]-b[i];w[i]=q[0]+(2*(p[i]-p[0])*(p[i]-p[0 ]))/(3*Pai);f=f+(Pai/60)*(q[i]-w[i])*(q[i]-w[i ])*(p[i]-p[i-1]);}return f;}void main(){floata,q,f,fl,f0,l[2],z[2],d0[100],d1[100], x[2],xi[2],fx,m0,m1,e;int i,j,n,k;printf("输入精度");scanf("%f",&e);do{z[0]=0+5*(rand()/;z[1]=0+5*(rand()/;}while(g(z[0],z[1])==0);for(i=0;i<=99;i++){d0[i]=-1+2*(rand()/;}for(j=0;j<=99;j++){d1[j]=-1+2*(rand()/;}f0=fun(z);fl=fun(z);ss:a=;for(i=0,j=0;i<=99&&j<=99;i++,j++) {n=1/sqrt((d0[i])*(d0[i])+d1[j]*d1[ j]);d0[i]=n*d0[i];d1[j]=n*d1[j];x[0]=z[0]+a*d0[i];x[1]=z[1]+a*d1[j];if(g(x[0],x[1])==1){f=fun(x);if(f<fl){fl=f;m0=d0[i];m1=d1[j];l[0]=x[0];l[1]=x[1];}}}x[0]=l[0];x[1]=l[1];do{a=*a;x[0]=x[0]+a*m0;x[1]=x[1]+a*m1;if(g(z[0],z[1])==0)break;f=fun(x);if(f<fl)fl=f;else break;} while(g(z[0],z[1])==1); do{x[0]=x[0]-a*m0;x[1]=x[1]-a*m1;a=*a;if(a<break;x[0]=x[0]+a*m0;x[1]=x[1]+a*m1; if(g(z[0],z[1])==1)f=fun(x);}while(f>=fl);if(fabs((f0-f)/f0)<e){xi[0]=x[0];xi[1]=x[1];fx=f;printf("最优解为\nx1*=%f\nx2*=%f\nfx=%f\n",xi[0],xi[1],fx);}else{f0=f;fl=f;z[0]=x[0];z[1]=x[1];goto ss;}}运行过程及结果：输入精度最优解为x1*=x2*=fx=。
机械优化设计外推法,黄金分割法,二次插值法
for(;a1>a3;) {t=a3; a3=a1; a1=t; t=y1; y3=y1;
if(f(x1)>f(x2)) *a=x1; else *b=x2; *n=*n+1; s=hj(a,b,e,n); } return s; } void main() { double s,a,b,e,m; int n=0; printf("输入 a，b 值和精度 e 值\n"); scanf("%lf %lf %lf",&a,&b,&e); s=hj(&a,&b,e,&n); m=(a+b)/2; printf("a=%lf,b=%lf,s=%lf,m=%lf,n=%d\n",a,b,s,m,n); }
} } while(sqrt(pow((x[2][1]-x[0][1]),2)+ pow((x[2][2]-x[0][2]),2))>=1e-6); xx[1]=x[2][1]; xx[2]=x[2][2]; fi=fun2(xx[1],xx[2]); printf("the best answer is : \nx1*= %f\nx2*=%f\nf*=%f\n",xx[1],xx[2],fi);
for(i=1;i<=2;i++) { if(i==1) {d[i][1]=1;
机械优化设计课件第三章-刘宇
厦门大学物理与机电工程学院
13
区间消去法原理
基本思想
搜索区间确定之后，采用区间消去法逐步缩短搜索区间，从而找 到极小点的数值近似解 在搜索区间[a, b]内任取两点a1和b1，且a1 < b1 ，计算其函数值得如 下结论 （1）f(a1 )<f(b1)，由于函数为单谷，所以极小点必在区间[a, b1]内 （2） f(a1 )>f(b1)，同理，极小点必在区间[a1, b]内 （3） f(a1 )=f(b1)，极小点在区间 极小点在区间[a1, b1]内
（1）选定初始点a1，初始步长h=h0，计算y1=f(α1)和y2=f(α1+h) （2）比较y1和y2 （a）如果y1>y2，转（3）向前探测 （b）如果y1<y2，h=-h0，将α1和α2，y1和y2的值互换。转（3）向后探测 （c）如果y1=y2，极小点在 极小点在α1和α1+h之间 （3）产生新的探测点α3= α2+h，y3=f(α3)
相同点：两种方法都是利用区间消去法原理将初始搜索区间不断 缩短，求得极小值的数值近似解 不同点：表现在试验点（插入点）位置的确定方法不同 试探法：试验点是按照某种个特定的规律确定；不考虑函数值的 试探法 试验点是按照某种个特定的规律确定 不考虑函数值的 分布 插值法：试验点是按照函数值近似分布的极小点确定；利用了函 数值本身及其导数信息；效果更好
第三章 一维搜索方法
主讲人：刘宇
2014/11/12
厦门大学物理与机电工程学院
1
本章内容
一、一维搜索的概念 维搜索的概念 二、搜索区间的确定与区间消去法原理 三 一维搜索的试探方法 三、 维搜索的试探方法—黄金分割法 四、一维搜索的插值方法
机械优化设计黄金分割法 外推法
郑州大学机械优化设计部分程序1、外推法2、黄金分割法3、二次插值法4、坐标轮换法5、随机方向法6、四杆机构优化设计1、外推法源程序：#include<stdio、h>#include<math、h>#define R 0、01double fun(double x){ double m;m=x*x-10*x+36;return m;}void main(){double h0=R,y1,y2,y3,x1,x2,x3,h; x1=0;h=h0;x2=h;y1=fun(x1);y2=fun(x2);if(y2>y1){h=-h;x3=x1;y3=y1;x1=x2;y1=y2;x2=x3;y2=y3;}x3=x2+h;y3=fun(x3); while(y3<y2){h*=2、0;x1=x2;y1=y2;x2=x3;y2=y3;x3=x2+h;y3=fun(x3);}printf("fun(%f)=%f,fun(%f)=%f,fun(%f)= %f\n",x1,y1,x2,y2,x3,y3);}运行过程及结果：fun(2、560000)=16、953600,fun(5、120000)=11、014400,fun(10、240000)=38、4576002、黄金分割法源程序：#include<stdio、h>#include<math、h>#define f(x) x*x*x*x-5*x*x*x+4*x*x-6*x+60double hj(double *a,double *b,double e,int *n) {double x1,x2,s;if(fabs((*b-*a)/(*b))<=e)s=f((*b+*a)/2);else{x1=*b-0、618*(*b-*a);x2=*a+0、618*(*b-*a);if(f(x1)>f(x2))*a=x1;else*b=x2;*n=*n+1;s=hj(a,b,e,n);}return s;}void main(){double s,a,b,e,m;int n=0;printf("输入a，b值与精度e值\n");scanf("%lf %lf %lf",&a,&b,&e);s=hj(&a,&b,e,&n);m=(a+b)/2;printf("a=%lf,b=%lf,s=%lf,m=%lf,n=%d\n",a,b ,s,m,n);}运行过程及结果：输入a，b值与精度e值-350、0001a=3、279466,b=3、279793,s=22、659008,m=3、279629,n=213、二次插值法源程序：#include<stdio、h>#include<math、h> int main(void){double a1,a2,a3,ap,y1,y2,y3,yp,c1,c2,m; double j[3];int i,h=1;void finding(double a[3]);finding(j);a1=j[0];a2=j[1];a3=j[2];m=0、001;double f(double x);y1=f(a1);y2=f(a2);y3=f(a3);for(i=1;1>=1;i++){c1=(y3-y1)/(a3-a1);c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);ap=0、5*(a1+a3-c1/c2);yp=f(ap);if(fabs((y2-yp)/y2)<m)break;else if((ap-a2)*h>0){if(y2>=yp){a1=a2;y1=y2;a2=ap;y2=yp;}else{a3=ap;y3=yp;}}else if(y2>=yp){a3=a2;y3=y2;a2=ap;y2=yp;}else{a1=ap;y1=yp;}}double x,y;if(y2<=yp){x=a2;y=y2;}else{x=ap;y=yp;}printf("a*=%f\n",x);printf("y*=%f\n",y);return 0;}double f(double x){double y;y=x*x-10*x+36;return y;}void finding(double a[3]){int h,i;double y[3];a[0]=0;h=1;a[1]=h;y[0]=f(a[0]);y[1]=f(a[1]);if(y[1]>y[0]){h=-h;a[2]=a[0];y[2]=y[0];do{a[0]=a[1];a[1]=a[2];y[0]=y[1];y[1]=y[2];a[2]=a[1]+h;y[2]=f(a[2]);h=2*h;}while(y[2]<y[1]);}else{for(i=1;i>=1;i++){a[2]=a[1]+h;y[2]=f(a[2]);if(y[2]>=y[1])break;h=2*h;a[0]=a[1];y[0]=y[1];a[1]=a[2];y[1]=y[2];}}return;}运行过程及结果：a*=5、000000y*=11、0000004、坐标轮换法源程序：#include <stdio、h>#include <math、h>#include <conio、h>float fun1(float x,float a,float b) {float y;y=x+a*b;return y;}float fun2(float x,float y){float z;z=4*(x-5)*(x-5)+(y-6)*(y-6);return z;}main(){float d[100][3],x[100][3],xx[3],ax[100][3]; float a1,a2,a3,h,t,y1,y2,y3,e,a,b,l,fi;int i,k;printf("输入初始点坐标\n");scanf("%f%f",&x[0][1],&x[0][2]);e=0、000001;l=0、618;x[2][1]=x[0][1];x[2][2]=x[0][2];k=0;k--;do{x[0][1]=x[2][1];x[0][2]=x[2][2];k++;for(i=1;i<=2;i++){if(i==1){d[i][1]=1;d[i][2]=0;}else{d[i][1]=0;d[i][2]=1;}h=0、1;a1=0;a2=h;x[i][1]=fun1(x[i-1][1],d[i][1],a1);x[i][2]=fun1(x[i-1][2],d[i][2],a1);y1=fun2(x[i][1],x[i][2]);x[i][1]=fun1(x[i-1][1],d[i][1],a2);x[i][2]=fun1(x[i-1][2],d[i][2],a2);y2=fun2(x[i][1],x[i][2]);if(y2>y1){h=-h;a3=a1;y3=y1;a1=a2;a2=a3;y1=y2;y2=y3;}a3=a2+h;x[i][1]=fun1(x[i-1][1],d[i][1],a3);x[i][2]=fun1(x[i-1][2],d[i][2],a3);y3=fun2(x[i][1],x[i][2]);do{a1=a2;y1=y2;a2=a3;y2=y3;a3=a2+h;x[i][1]=fun1(x[i-1][1],d[i][1],a3);x[i][2]=fun1(x[i-1][2],d[i][2],a3);y3=fun2(x[i][1],x[i][2]);}while(y3<y2);for(;a1>a3;){t=a3;a3=a1;a1=t;t=y1;y3=y1;y1=t;}a=a1;b=a3;a1=b-l*(b-a);a2=a+l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a1);x[i][2]=fun1(x[i-1][2],d[i][2],a1);y1=fun2(x[i][1],x[i][2]);x[i][1]=fun1(x[i-1][1],d[i][1],a2);x[i][2]=fun1(x[i-1][2],d[i][2],a2);y2=fun2(x[i][1],x[i][2]);if(b<1e-3){for(;fabs(b-a)>e;){if(y1>=y2){a=a1;a1=a2;y1=y2;a2=a+l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a2);x[i][2]=fun1(x[i-1][2],d[i][2],a2);y2=fun2(x[i][1],x[i][2]);}else{b=a2;a2=a1;y2=y1;a1=b-l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a1);x[i][2]=fun1(x[i-1][2],d[i][2],a1);y1=fun2(x[i][1],x[i][2]);}}}else{for(;fabs((b-a)/b)>=e||fabs((y2-y1)/y2)>=e; ){if(y1>=y2){a=a1;a1=a2;y1=y2;a2=a+l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a2);x[i][2]=fun1(x[i-1][2],d[i][2],a2);y2=fun2(x[i][1],x[i][2]);}else{b=a2;a2=a1;y2=y1;a1=b-l*(b-a);x[i][1]=fun1(x[i-1][1],d[i][1],a1);x[i][2]=fun1(x[i-1][2],d[i][2],a1);y1=fun2(x[i][1],x[i][2]);}}}ax[k][i]=0、5*(a+b);x[i][1]=fun1(x[i-1][1],d[i][1],ax[k][i]);x[i][2]=fun1(x[i-1][2],d[i][2],ax[k][i]);}}while(sqrt(pow((x[2][1]-x[0][1]),2)+pow((x[2][2]-x[0][2]),2))>=1e-6);xx[1]=x[2][1];xx[2]=x[2][2];fi=fun2(xx[1],xx[2]);printf("最优解为\nx1*=%f\nx2*=%f\nf*=%f\nk=%d\n",xx[1],xx[2],fi,k);}运行过程及结果：输入初始点坐标89最优解为x1*=5、000000x2*=6、000000f*=0、000000k=25、随机方向法源程序：#include<math、h>#include<stdio、h>#include<stdlib、h>float f(float x,float y){float z;z=(x-2)*(x-2)+(y-1)*(y-1);return z;}float g1(float x,float y){float z;z=x*x-y;return z;}float g2(float x,float y){float z;z=x+y-2;return z; }void main(){int i,j;float k=8,c=0、000001,a0=-3,b0=3,a1=-3,b1=3;floatx[10],x0[10],xl[10],e[10],r[10],d[10],h,fl,f0,fx; while(g1(x0[0],x0[1])>0||g2(x0[0],x0[1])>0) {x0[0]=a0+(rand()/32767、00)*(b0-a0);x0[1]=a1+(rand()/32767、00)*(b1-a1);}fl=f(x0[0],x0[1]);f0=f(x0[0],x0[1]);while(1){h=0、01;j=1;r[0]=-1+(rand()/32767、00)*(1-(-1));r[1]=-1+(rand()/32767、00)*(1-(-1));e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);x[0]=x0[0]+h*e[0];x[1]=x0[1]+h*e[1];if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0){fx=f(x[0],x[1]);if(fx<fl){fl=fx;for(i=0;i<2;i++){d[i]=e[i];xl[i]=x[i];}}}while(j<=k){j++;r[0]=-1+(rand()/32767、00)*(1-(-1));r[1]=-1+(rand()/32767、00)*(1-(-1));e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);x[0]=x0[0]+h*e[0];x[1]=x0[1]+h*e[1];if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0){fx=f(x[0],x[1]);if(fx<fl){fl=fx;for(i=0;i<2;i++){d[i]=e[i];xl[i]=x[i];}}}}x[0]=xl[0];x[1]=xl[1];while(1){h=1、3*h;x[0]=x[0]+h*d[0];x[1]=x[1]+h*d[1];if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)break;fx=f(x[0],x[1]);if(fx<fl) fl=fx;else break;}do{x[0]=x[0]-h*d[0];x[1]=x[1]-h*d[1];h=0、7*h;if(h<c)break;x[0]=x[0]+h*d[0];x[1]=x[1]+h*d[1];if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)continue;fx=f(x[0],x[1]);}while(fx>=fl);if(fabs((f0-fx)/f0)>=c){x0[0]=x[0];x0[1]=x[1];fl=fx;f0=fx;}elsebreak;}printf("输出最优解为\nx1*=%f,x2*=%f,y*=%f\n",x[0],x[1],fx);}运行过程及结果：输出最优解为x1*=0、995421,x2*=1、004521,y*=1、009200 6、四杆机构优化设计源程序：#include<math、h>#include<stdio、h>#include<stdlib、h>#define Pai 3、1415926int g(float l1,float l2){if((-l1<=0)&&(-l2<=0)&&(6-l1-l2<=0)&&(1-l2-4<=0)&&(l2-l1-4<=0)&&(l1*l1+l2*l2-1、414*l1*l2-16<=0)&&(36-l1*l1-l2*l2-1、414*l1*l2<=0))return (1);elsereturn (0);}float fun(float x0[2]){floatf,a[31],b[31],r[31],p[31],q[31],w[31],x1[2];int i;p[0]=acos(((1+x0[0])*(1+x0[0])-x0[1]*x0[ 1]+25)/(10+10*x0[0]));q[0]=acos(((1+x0[0])*(1+x0[0])-x0[1]*x0[ 1]-25)/(10*x0[1]));f=0;for(i=1;i<=30;i++){p[i]=p[0]+(Pai/60)*i;r[i]=sqrt(26-10*cos(p[i]));a[i]=acos((r[i]*r[i]+x0[1]*x0[1]-x0[0]*x0[ 0])/(2*r[i]*x0[1]));b[i]=acos((r[i]*r[i]+24)/(10*r[i]));q[i]=Pai-a[i]-b[i];w[i]=q[0]+(2*(p[i]-p[0])*(p[i]-p[0]))/(3*P ai);f=f+(Pai/60)*(q[i]-w[i])*(q[i]-w[i])*(p[i]-p [i-1]);}return f;}void main(){floata,q,f,fl,f0,l[2],z[2],d0[100],d1[100],x[2],xi[2],f x,m0,m1,e;int i,j,n,k;printf("输入精度");scanf("%f",&e);do{z[0]=0+5*(rand()/32767、00);z[1]=0+5*(rand()/32767、00);}while(g(z[0],z[1])==0);for(i=0;i<=99;i++){d0[i]=-1+2*(rand()/32767、00);}for(j=0;j<=99;j++){d1[j]=-1+2*(rand()/32767、00);}f0=fun(z);fl=fun(z);ss:a=0、01;for(i=0,j=0;i<=99&&j<=99;i++,j++) {n=1/sqrt((d0[i])*(d0[i])+d1[j]*d1[j]);d0[i]=n*d0[i];d1[j]=n*d1[j];x[0]=z[0]+a*d0[i];x[1]=z[1]+a*d1[j];if(g(x[0],x[1])==1){f=fun(x);if(f<fl){fl=f;m0=d0[i];m1=d1[j];l[0]=x[0];l[1]=x[1];}}}x[0]=l[0];x[1]=l[1];do{a=1、3*a;x[0]=x[0]+a*m0;x[1]=x[1]+a*m1;if(g(z[0],z[1])==0)break;f=fun(x);if(f<fl)fl=f;else break;}while(g(z[0],z[1])==1);do{x[0]=x[0]-a*m0;x[1]=x[1]-a*m1;a=0、7*a;if(a<0、00001)break;x[0]=x[0]+a*m0;x[1]=x[1]+a*m1;if(g(z[0],z[1])==1)f=fun(x);}while(f>=fl);if(fabs((f0-f)/f0)<e){xi[0]=x[0];xi[1]=x[1];fx=f;printf("最优解为\nx1*=%f\nx2*=%f\n fx=%f\n",xi[0],xi[1],fx);}else{f0=f;fl=f;z[0]=x[0];z[1]=x[1];goto ss;}}运行过程及结果：输入精度0、001最优解为x1*=4、161386x2*=2、311257fx=0、000021。
机械优化设计-第02章 优化计算方法
第02章优化计算方法2.1黄金分割法黄金分割法也称0.618法，是通过对黄金分割点函数值的计算和比较，将初始区间逐次进行缩小，直到满足给定的精度要求，即求得一维极小点的近似解。
一、方法概述（一）区间缩小的基本思路已知的单峰区间。
为了缩小区间，在内按一定规则对称地取2个内部点和，并计算和。
可能有三种情况：图(a)经过一次函数比较，区间缩小一次。
在新的区间内，保留一个好点和，下一次只需再按一定规则，在新区间内找另一个与对称的点，计算，与比较。
如此反复。
图(b)淘汰，另，得新区间。
图(c)可归纳入上面任一种情况处理。
（二）取点规则黄金分割法的关键是如何不断找出区间内的2个对称点，保证极小点不会丢掉，且收敛快。
设初始区间长度为l，第一次区间缩短率为，则缩短后的区间长度为。
第二次区间缩短时，在区间中取点，经比较后又得新区间。
由对称性可知，区间的长度为，则本次区间缩短率为令这两次缩短率相等，即，得方程解方程，得合理的根为由此可知，黄金分割法的均匀缩短率为0.618，即每经过一次函数值比较，都是淘汰本次区间的0.382倍。
根据上式，黄金分割法的取点规则是为了使最终区间收敛到给定收敛精度内，区间的缩短次数N必须满足：即二、收敛准则由于实际问题的需要和函数形态的不同，常常需要不同的收敛准则确定最优点。
对于直接法，有以下几种收敛准则：（1）区间绝对精度；（2）区间相对精度；（3）函数值绝对精度；（4）函数值相对精度三、方法特点（一）黄金分割法特点（1）不必要求可微，只要利用函数值大小的比较，即可很快地找到；（2）除了第一次缩小区间要计算两个点及其函数值以外，其余每次只要计算一个点及其函数值；（3）可靠性好。
（二）应用举例实际一个圆柱螺旋压缩弹簧，不考虑共振，要求重量W最轻。
解：建模前，先列出弹簧的有关设计计算公式：式中-------弹簧的设计载荷；-------弹簧的总变形量；-------弹簧指数；K-------曲度系数；n-------工作有效圈数；n2-------不起作用圈数（总圈数与工作有效圈数之差）；-------材料密度。
4 现代设计方法--优化设计
ADM
目录
第三章 平面问题有限元 3.1 平面问题基本方程及有限元矩阵方程 3.1.1 基本方程 3.1.2 有限元矩阵方程 3.2 三角形场应变单元 3.2.1 离散化 3.2.2 位移模式 3.2.3 应变 3.4 刚度矩阵 3.4.1 单元刚度矩阵 3.4.2 总体刚度矩阵的组装 3.4.3 总体位移向量 3.5 单元的等效节点力与总体载荷向量 3.5.1 单元的等效节点力 3.5.2 总体载荷向量
现代设计方法
——优化设计、有限元
Advanced Design Methods
——Design Optimization and Finite Element Method
江南大学 机械工程学院
1
ADM
目录
序论
第一部分 优化设计
第一章 优化设计的数学基础
1.1 矢量 1.2 矩阵 1.3 多元函数
目录
7
ADM
目录
第六章 杆件系统 第七章 薄板弯曲问题 第八章 结构动力学问题
8.1 结构动力学微分方程 8.2 结构动力学虚功方程 8.3 结构动力学有限元矩阵方程 8.4 结构自由振动有限元矩阵方程——模态分析
8
ADM
序论
现代设计方法的基本内容：
1. CAD 2. CAE——有限元分析* 3. 优化设计* 4. 可靠性设计 5. 逆向设计 6. 模块化设计 7. 设计专家系统 8. 价值工程 9. 虚拟设计 10. ……………
F(X0) 0
极值存在的充分条件：
DF
DX TF(X0 )
1 2
DX T H (X0 )DX
1 2
DX T H (X0 )DX
H(X0)正定， F(X0)为极小值；