机械优化设计 变尺度法 c语言程序

计算 f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2 的无约束极值，唤掘初始点x0=[1,1]。

/*

tt ---- 一维搜索初始步长

ff ---- 差分法求梯度时的步长岩链迟

ac ---- 终止迭代收敛精度

ad ---- 一维搜索收敛精度

n ----- 设计变量的维数

xk[n] -- 迭代初始点

*/

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#include<conio.h>

#define tt 0.01

#define ff 1.0e-6

#define ac 1.0e-6

#define ad 1.0e-6

#define n 2

double ia

double fny(double *x)

{

double x1=x[0],x2=x[1]

double f

f=x1*x1+2*x2*x2-4*x1-2*x1*x2

return f

}

double * iterate(double *x,double a,double *s)

{

double *x1

int i

x1=(double *)malloc(n*sizeof(double))

for(i=0i<ni++)

x1[i]=x[i]+a*s[i]

return x1

}

double func(double *x,double a,double *s)

{

double *x1

double f

x1=iterate(x,a,s)

f=fny(x1)

return f

}

void finding(double a[3],double f[3],double *xk,double *s)

{

double t=tt

int i

double a1,f1

a[0]=0f[0]=func(xk,a[0],s)

for(i=0i++)

{

a[1]=a[0]+t

f[1]=func(xk,a[1],s)

if(f[1]<f[0]) break

if(fabs(f[1]-f[0])>=ad)

{

t=-t

a[0]=a[1]f[0]=f[1]

}

else

{

if(ia==1) return//break

t=t/2ia=1

}

}

for(i=0i++)

{

a[2]=a[1]+t

f[2]=func(xk,a[2],s)

if(f[2]>f[1]) break

t=2*t

a[0]=a[1]f[0]=f[1]

a[1]=a[2]f[1]=f[2]

}

if(a[0]>a[2])

{

a1=a[0]

f1=f[0]

a[0]=a[2]

f[0]=f[2]

a[2]=a1

f[2]=f1

}

return

}

double lagrange(double *xk,double *ft,double *s)

{

int i

double a[3],f[3]

double b,c,d,aa

finding(a,f,xk,s)

for(i=0i++)

{

if(ia==1) { aa=a[1]*ft=f[1]break}

d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2])

if(fabs(d)==0) break

c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d

if(fabs(c)==0) break

b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1])

aa=-b/(2*c)

*ft=func(xk,aa,s)

if(fabs(aa-a[1])<=ad) {if(*ft>f[1]) aa=a[1]break}

if(aa>a[1])

{

if(*ft>粗李f[1]) {a[2]=aaf[2]=*ft}

else if(*ft<f[1]) {a[0]=a[1]a[1]=aaf[0]=f[1]f[1]=*ft}

else if(*ft==f[1])

{

a[2]=aaa[0]=a[1]

f[2]=*ftf[0]=f[1]

a[1]=(a[0]+a[2])/2

f[1]=func(xk,a[1],s)

}

}

else

{

if(*ft>f[1]) {a[0]=aaf[0]=*ft}

else if(*ft<f[1]) {a[2]=a[1]a[1]=aaf[2]=f[1]f[1]=*ft}

else if(*ft==f[1])

{a[0]=aaa[2]=a[1]

f[0]=*ftf[2]=f[1]

a[1]=(a[0]+a[2])/2

f[1]=func(xk,a[1],s)

}

}

}

if(*ft>f[1]) {*ft=f[1]aa=a[1]}

return aa

}

double *gradient(double *xk)

{

double *g,f1,f2,q

int i

g=(double*)malloc(n*sizeof(double))

f1=fny(xk)

for(i=0i<ni++)

{q=ff

xk[i]=xk[i]+qf2=fny(xk)

g[i]=(f2-f1)/qxk[i]=xk[i]-q

}

return g

}

double * bfgs(double *xk)

{

double u[n],v[n],h[n][n],dx[n],dg[n],s[n]

double aa,ib

double *ft,*xk1,*g1,*g2,*xx,*x0=xk

double fi

int i,j,k

ft=(double *)malloc(sizeof(double))

xk1=(double *)malloc(n*sizeof(double))

for(i=0i<ni++)

{

s[i]=0

for(j=0j<nj++)

{

h[i][j]=0

if(j==i) h[i][j]=1

}

}

g1=gradient(xk)

fi=fny(xk)

x0=xk

for(k=0k<nk++)

{

ib=0

if(ia==1) { xx=xkbreak}

ib=0

for(i=0i<ni++) s[i]=0

for(i=0i<ni++)

for(j=0j<nj++)

s[i]+= -h[i][j]*g1[j]

aa=lagrange(xk,ft,s)

xk1=iterate(xk,aa,s)

g2=gradient(xk1)

for(i=0i<ni++)

if((fabs(g2[i])>=ac)&&(fabs(g2[i]-g1[i])>=ac))

{ib=ib+1}

if(ib==0) { xx=xk1break}

fi=*ft

if(k==n-1)

{ int j

xk=xk1

for(i=0i<ni++)

for(j=0j<nj++)

{

h[i][j]=0

if(j==i) h[i][j]=1

}

g1=g2k=-1

}

else

{

int j

double a1=0,a2=0

for(i=0i<ni++)

{

dg[i]=g2[i]-g1[i]

dx[i]=xk1[i]-xk[i]

}

for(i=0i<ni++)

{

int j

u[i]=0v[i]=0

for(j=0j<nj++)

{

u[i]=u[i]+dg[j]*h[j][i]

v[i]=v[i]+dg[j]*h[i][j]

}

}

for(j=0j<nj++)

{

a1+=dx[j]*dg[j]

a2+=v[j]*dg[j]

}

if(fabs(a1)!=0)

{

a2=1+a2/a1

for(i=0i<ni++)

for(j=0j<nj++)

h[i][j]+=(a2*dx[i]*dx[j]-v[i]*dx[j]-dx[i]*u[j])/a1

}

xk=xk1g1=g2

}

}

if(*ft>fi) { *ft=fixx=xk}

xk=x0

return xx

}

void main ()

{

int k

double *xx,f

double xk[n]={1,1}

xx=bfgs(xk)

f=fny(xx)

printf("\n\nThe Optimal Design Result Is:\n")

for(k=0k<nk++)

{printf("\n\tx[%d]*=%f",k+1,xx[k])}

printf("\n\tf*=%f",f)

getch()

}

这是基于一本书上的算法。但我很奇怪，原书中的算法有结果列出，但是我却不能通过编译。真是纳闷！修改后可以得到结果了，如果你要使用这个简单的程序，你只需更改 维数n、double fny(double *x)的实现部分以及main函数中的xk初值就可以了。不过这个程序也不是很好。

>>clear

>>f=inline('a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4)','a','x')

x=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]

y=[2.3201 2.6470 2.9707 3.2885 3.6008 3.9090 4.2147 4.5191 4.8232 5.1275]

[xx,res]=lsqcurvefit(f,ones(1,4),x,y)

xx'宽梁,res

要建立也是可以的。就是空裂把上面那个inline弄成.m

如下：

在Matlab下输入：edit zhidao_15.m，然后将下面两行百分号之间的斗巧闭内容，复制进去，保存

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y=zhidao_15(a,x)

y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

在Matlab下面输入：

x=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]

y=[2.3201 2.6470 2.9707 3.2885 3.6008 3.9090 4.2147 4.5191 4.8232 5.1275]

[xx,res]=lsqcurvefit('zhidao_15',ones(1,4),x,y)

xx',res

一、概观

scipy中的optimize子包中提供了常用的最优化算法函数实现。我们可以直接调用这些函数完成我们的优化问题。optimize中函数最典型的特点就是能够从函数名称上看出是使用了什么算法。下面optimize包中函数的概览：

1.非线性最优化

fmin -- 简单Nelder-Mead算法悔枯

fmin_powell -- 改进型Powell法

fmin_bfgs -- 拟Newton法

fmin_cg -- 非线性共轭梯度法

fmin_ncg -- 线性搜索Newton共轭梯度法

leastsq -- 最小二乘

2.有约束的多元函数问题

fmin_l_bfgs_b ---使用L-BFGS-B算法

fmin_tnc ---梯度信息

fmin_cobyla ---线性逼近

fmin_slsqp ---序列最小二乘法

nnls ---解|| Ax - b ||_2 for x>=0

3.全局优化

anneal ---模拟退火算法

brute --强力法

4.标量函数

fminbound

brent

golden

bracket

5.拟合

curve_fit-- 使用非线性最小二乘法拟合

6.标量函数求根

brentq ---classic Brent (1973)

brenth ---A variation on the classic Brent（1980）ridder ---Ridder是提出这个算法的人培搏名

bisect ---二分法

newton ---牛顿法

fixed_point

7.多维函数求根

fsolve ---通用

broyden1 ---Broyden’s first Jacobian approximation.

broyden2 ---Broyden’s second Jacobian approximationnewton_krylov ---Krylov approximation for inverse Jacobiananderson ---extended Anderson mixing

excitingmixing ---tuned diagonal Jacobian approximationlinearmixing ---scalar Jacobian approximationdiagbroyden ---diagonal Broyden Jacobian approximation8.实用函数

line_search ---找到满足强Wolfe的alpha值

check_grad ---通过和前向有限差分逼近比较检查梯度函数的正确性二、实战非线性最优化

fmin完整的调用形式是：

fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None)不过我们最常使用的就是前两个参数。一个描述优化问题的函数以及初值。后面的那些参数我们也很容易理解。如果您能用到，请自己研究。下面研究一个最简单的问题，来感受这个函数的使用方法：f(x)=x**2-4*x+8，我们知道，这个函数的最小值是4，在x=2的时候取到。

from scipy.optimize import fmin#引入优化包def myfunc(x):

return x**2-4*x+8#定义函数

x0 = [1.3]#猜一个初值

xopt = fmin(myfunc, x0)#求解

print xopt#打印结果

运行之后，给出的结果是：

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 16

Function evaluations: 32

[ 2.00001953]

程序准确的计算得出了最小值，不过最小值点并不是严格的2，碧中洞这应该是由二进制机器编码误差造成的。

除了fmin_ncg必须提供梯度信息外，其他几个函数的调用大同小异，完全类似。我们不妨做一个对比：

from scipy.optimize import fmin,fmin_powell,fmin_bfgs,fmin_cgdef myfunc(x):

return x**2-4*x+8

x0 = [1.3]

xopt1 = fmin(myfunc, x0)

print xopt1

print

xopt2 = fmin_powell(myfunc, x0)

print xopt2

print

xopt3 = fmin_bfgs(myfunc, x0)

print xopt3

print

xopt4 = fmin_cg(myfunc,x0)

print xopt4

给出的结果是：

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 16

Function evaluations: 32

[ 2.00001953]

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 2

Function evaluations: 53

1.99999999997

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 2

Function evaluations: 12

Gradient evaluations: 4

[ 2.00000001]

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 2

Function evaluations: 15

Gradient evaluations: 5

[ 2.]

我们可以根据给出的消息直观的判断算法的执行情况。每一种算法数学上的问题，请自己看书学习。个人感觉，如果不是纯研究数学的工作，没必要搞清楚那些推导以及定理云云。不过，必须了解每一种算法的优劣以及能力所及。在使用的时候，不妨多种算法都使用一下，看看效果分别如何，同时，还可以互相印证算法失效的问题。

在from scipy.optimize import fmin之后，就可以使用help(fmin)来查看fmin的帮助信息了。帮助信息中没有例子，但是给出了每一个参数的含义说明，这是调用函数时候的最有价值参考。

有源码研究癖好的，或者当你需要改进这些已经实现的算法的时候，可能需要查看optimize中的每种算法的源代码。在这里：https:/ / github. com/scipy/scipy/blob/master/scipy/optimize/optimize.py聪明的你肯定发现了，顺着这个链接往上一级、再往上一级，你会找到scipy的几乎所有源码！

---