最优化大作业2014

题目：最优化大作业班级：021251

姓名：杨明坤

学号：02125043

：无论做什么事，人们总希望以最少的代价取得最大的利益，也就是力求最好，这就是最优化问题。最优化就是在一切可能的方案中，选择一个最好的方案以达到最优目标。

概括地说，凡是追求最有目标的数学问题都属于最优化问题。作为最优化问题，一般有三个要素：

第一是目标

第二是方案

第三是限制条件

而且目标应是方案的“函数”，如果方案与时间无关，则该问题属于静态最优化问题;否则属于动态最优化问题。

爬坡法求解y=(r-2）*（r-2)+3的局部最优解1.首先介绍一些基本概念：

上面是爬山法的图示描述。下面给出本体的实现代码: #include

#include

int com(float,float);

float x,i,z,w,h;

printf("%f",&x);

main()

{

z=x=0;

i=0.1;

for(i<=6)

{

w=x+i;

h=com(x,w);

if(h==1);

z=w;

x=w;

}

printf("%f/n",f(w));

}

float f(float r)

{

float y;

y=(r-2)*(r-2)+3;

return y;

}

int com(float x ,float i)

{floata,b;

a=f(x);

b=f(i);

if((b-a)<0);

return 1;

else

return o;

}

下面是最优化问题的数学模型及分类;

用修正牛顿法求解minf(X)=4(x1+1)^2+2(x2-1)^2+x1+x2+10,初始点X0=【0,0】^T,?=0.01. 先介绍一些基本概念：

正定二次函数，牛顿反向就是指极小点的方向。

用牛顿法解目标函数的无约束问题，只需一步迭代就可以得到最优解。

下面介绍修正牛顿法:

说明：修正牛顿法对初始点的选取要求不严。

下面是解题步骤;

源代码：

#include

#include

#include

#define n 2//正定二次函数的自变量个数

double fun(double x[n],double f_xs[n+n+1+(n-1)*n/2])//输入变量为函数自变量初值{

int i,j;

double f=0;

for(i=0;i

{

f+=pow(x[i],2)*f_xs[i];

}

for(;i<2*n;i++)//计算X部分

{

f+=x[i%n]*f_xs[i];

}

f+=f_xs[i];//计算常数项部分

for(i=0;i

{

for(j=i+1;j

{

f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];

}

}

return f;

}

void Q_fun(double f_xs[n+n+1+(n-1)*n/2],double Q[n][n+1])

{

int i,j;

for(i=0;i

{

Q[i][i]=2*f_xs[i];

}

for(i=0;i

{

Q[i][n]=f_xs[n+i];

}

for(i=0;i

{

for(j=i+1;j

{

Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];

}

}

}

void D_fun(double x[n],double Q[n][n+1],double g0[n])//自变量初值，正定二次函数的各项系数，返回梯度的初值

{

int i,j;

for(i=0;i

{

g0[i]=0;//清零

for(j=0;j

{

if(i==j)

g0[i]+=Q[i][j]*x[j];

else

g0[i]+=Q[i][j]*x[j];

}

}

for(i=0;i

{

g0[i]+=Q[i][n];

}

cout<

for(i=0;i

cout<

cout<

}

void Hesse_fun(double Q[n][n+1],double Hesse[n][n])

{

int i,j;

for(i=0;i

{

for(j=0;j

{

Hesse[i][j]=Q[i][j];

}

}

cout<<"Hesse 矩阵:"<

for(i=0;i

{

for(j=0;j

{

cout<

}

cout<

}

}

int H(double g0[n],double c)//判别准则：返回1结束，返回0继续迭代{

double s=0;

for(int i=0;i

{

s+=pow(g0[i],2);

}

if(sqrt(s)

return 1;

else

return 0;

}

void inv(double a[n][n],double c[n][n])//求Hesse矩阵的逆矩阵

{

double m[n];//辅助乘数

double T;//存储换行临时变量

double temp=0;

double t;//最大列主元

int tap1=0,tap2=0;//最大列主元下标

for(int i=0;i

{

for(int j=0;j

{

if(i==j)

{

c[i][j]=1;

}

else

{

c[i][j]=0;

}

}

}

for(int s=0;s

{

t=a[s][s];//赋初值

for(int p=s;p

{

if(fabs(a[p][s])>t)

{

t=a[p][s];

tap1=p;

tap2=s;

}

}

if(t==0)

{cout<<"列主元为零，失败！"<

{

for(int i=0;i

{

T=a[tap1][i];

a[tap1][i]=a[s][i];

a[s][i]=T;

T=c[tap1][i];//逆矩阵换行

c[tap1][i]=c[s][i];

c[s][i]=T;

}

}

tap1=tap2=0;//置零

for(int j=0;j

{

a[s][j]=a[s][j]/t;

c[s][j]=c[s][j]/t;

}

for(int i=0;i

{

if(i!=s)

{

m[i]=a[i][s]/a[s][s];

for(int j=0;j

{

a[i][j]=a[i][j]-m[i]*a[s][j];

c[i][j]=c[i][j]-m[i]*c[s][j];

}

}

}

}

cout<

cout<<"Hesse矩阵的逆矩阵为:"<

for(i=0;i

{

for(int j=0;j

{

cout<

}

cout<

}

}

void xiang_cheng(double a[n][n],double b[n],double c[n])//n*n矩阵乘以n*1矩阵

{

int i,j,k;

if(n==n)

{

cout<<"两矩阵阶数相符可以相乘！"<

for(i=0;i

{

c[i]=0;

}

for(k=0;k

{

for(j=0;j

{

c[k]+=a[k][j]*b[j];

}

}

cout<<"Hesse矩阵的逆矩阵和梯度向量相乘，两矩阵相乘结果为:"<

for(i=0;i

{

cout<

}

}

else

cout<<"两个矩阵阶数不符，不能相乘"<

}

//开始计算minf(Xk+tPk)时的步长t的值，由于这是n元二次函数所以minf(t)是关于t>0的二次函数，先求二次方程a,b,c系数,用一阶导为零。

void abc(double x[n],double p[n],double f_xs[n+n+1+(n-1)*n/2],double t[3])//t[3]中返回的是a,b,c的系数值

{

int i,j;

t[0]=t[1]=t[2]=0;

t[0]=fun(p,f_xs)-f_xs[2*n];

for(i=n;i<2*n;i++)

{

t[0]-=f_xs[i]*p[i%n];

}

for(i=0;i

{

t[1]+=2*f_xs[i]*x[i]*p[i];

}

for(;i<2*n;i++)

{

t[1]+=f_xs[i]*p[i%n];

}

for(i=0;i

{

for(j=i+1;j

{

t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]);

}

}

t[2]=fun(x,f_xs);

cout<

}

//二次函数一阶导为零计算t的值，t>0

double t_c(double t[3])

{

cout<<"用一阶导求得步长为：t="<<-t[1]/t[0]/2<

return -t[1]/t[0]/2;

}

void main()

{

double f_xs[n+n+1+(n-1)*n/2]={4,2,9,-3,16,0};//正定二次函数的各项系数，第一个为：X1^2系数，第二个为：X2^2系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数；

//n元的系数存放类推

double x[n]={0,0};//函数自变量初值

double f0;//函数值

double g0[n];//梯度的值

double Q[n][n+1];//正定二次函数对应的实对称矩阵

double c=0.01;//H准则限值

double Hesse[n][n];//Hesse矩阵

double inv_Hesse[n][n];//Hesse矩阵的逆

double t[3];//返回求minf()时t的二次函数的a,b,c的系数值

double t_bc;//步长

double p[n];//保存矩阵相乘结果_加负号后变成下降方向

int i,tap=0;//迭代次数

Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵ok!

f0=fun(x,f_xs);//求函数初值ok!

D_fun(x,Q,g0);//返回梯度的初值

while(H(g0,c)==0)

{

//x_k_1(x,g0,Q);//求出下一个迭代点，更新了x[n]的值

Hesse_fun(Q,Hesse);//返回Hesse矩阵的值？？一个二次函数的Hesse矩阵是变的还是不变

inv(Hesse,inv_Hesse);//返回Hesse矩阵的逆矩阵

xiang_cheng(inv_Hesse,g0,p);

for(i=0;i

p[i]*=(-1);

abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值，

t_bc=t_c(t);//求一阶导来计算t

for(i=0;i

{

x[i]=x[i]+t_bc*p[i];

}

tap++;

f0=fun(x,f_xs);

D_fun(x,Q,g0);

}

cout<<"修正newton法"<

cout<<"函数f(x1,x2)=4(X1+1)^2+2(X2-1)^2+X1+X2+10.的极小点为："<<"f="<

cout<<"自变量取值为："<

for(i=0;i

{

cout<<"x"<

}

cout<<"迭代次数为："<

}

用外点罚函数法编程计算;

Minf(X)=-x1+x2

S.t.g1(X)=lnx2≥0

h1(X)=x1+x2-1=0

取终止限ε=10^-5

下面先介绍基本知识：

外点法;

下面是迭代过程;

最后附上源代码;

function F=fun(x1,x2,m,u)

f=-x1+x2;

g=x2-1;

h=x1+x2-1;

F=f+m*h.^2+m*g.^2*u;

clear all

clc

m=1;c=10;d=10^(-15);

k=1;

fw(1)=fun(3,4,0,0);xw(:,1)=[3;4];

while 1

x(:,1)=[3;4];

if (x(2,1)-1>=0)

u=0;

else

u=1;

end

syms x1x2

ff=fun(x1,x2,m,u);

gg=jacobian(ff,[x1;x2])';

GG=jacobian([gg(1);gg(2)],[x1;x2]);

%=====================================

%一些初始值的计算

f(1)=subs(ff,{x1,x2},{x(1,1),x(2,1)});

g=subs(gg,{x1,x2},{x(1,1),x(2,1)});

G=subs(GG,{x1,x2},{x(1,1),x(2,1)});

i=1;

e=10;

while e>=10^(-10)

t(:,i)=(g'*g/(g'*G*g))*g;

x(:,i+1)=x(:,i)-t(:,i);

f(i+1)=fun(x(1,i+1),x(2,i+1),m,u);

g=subs(gg,{x1,x2},{x(1,i+1),x(2,i+1)});

G=subs(GG,{x1,x2},{x(1,i+1),x(2,i+1)});

e=f(i+1)-f(i);

i=i+1;

end

%=================================================== k=k+1;

fw(k)=f(i);xw(:,k)=x(:,i);

%//////////////////////////////////////////////////////////////// ////////////////

if abs(m*(x(1,i)+x(2,i)-1).^2+m*u*(x(2,i)-1).^2)<=d

disp('迭代过程X');

xw

fw

break;

else

m=m*c;

end

end