实验的目的和要求：通过本次实验使学生较为熟练使用MATLAB软件，并能利用该软件进行约束最优化方法的计算。
实验内容：
１、罚函数法的MATLAB实现
２、可行方向法的MATLAB实现
学习建议：
本次实验就是要通过对一些具体问题的分析进一步熟悉软件的操作并加深对理论知识的理解。
重点和难点：
可行点和辅助函数选取。
一 罚函数法
用罚函数法解min f(x)=(x1-1)2+x22
S.T. g(x)=x2-1>=0
编写下面m文件
fahanshu.m
syms x1 x2
f=(x1-1)^2+x2^2;
g=x2-1;
x=[x1;x2];
x0=[0;0];
e=0.0001;M=1;
while abs(subs(g,x,x0))>e
if subs(g,x,x0)<0
Q=f+M*g^2;
else
Q=f;
end
x0=xzNewton(Q,x,x0,0.0001);
M=10*M;
end
x0
运行得：fahanshu
x0 =
1.0000
0.9999
与理论值x=[1;1]很接近。
二 投影梯度法
1.梯度投影法基本原理和步骤
思想：当迭代点是可行域的内点时，将目标函数负梯度作为搜索方向，当迭代点在可行域边界上时，将目标函数负梯度在可行域边界上的投影作为搜索方向。无论何种情况，所构造的方向都是可行下降方向。然后在可行域内沿该方向进行最优一维搜索得到新的迭代点。
MATLAB实现:
2.代码及数值算例:
(1) 程序源代码:
function [ X,FMIN,K ] = tidutouying( f,A,b,x1,x,e )
% [ X,FMIN,K ] = tidutouying( f,A,b,x1,e ) 梯度投影法
% f 目标函数
% A 约束矩阵 b 右端项
% x1 初始点 x 自由变量
% e 精度要求
% X 极小点
% FMIN 极小值
% K 迭代次数
% 张超编写与2014/5/3
count=1;
n=length(x1);
tf=jacobian(f,x)';
while 1
[A1,A2,b1,b2,k]=fenjie(A,b,x1);
while 1
M=A1;
if isempty(M)
P=eye(n);
else
P=eye(n)-M'*(M*M')^(-1)*M;
end
Pk=-P*subs(tf,x,x1);
if norm(Pk)<=e
if isempty(M)
x1;break;
else
W=(M*M')^(-1)*M*subs(tf,x,x1);
u=W;
if min(u)<0
for i=1:length(u)
if u(i)==min(u)
j=i;
end
end
A1(j,:)=[];
else
x1;break;
end
end
else
b_=b2-A2*x1;
P_=A2*Pk;
for i=1:length(P_)
if P_(i)<0
r(i)=b_(i)/P_(i);
else
r(i)=10000;
end
end
rmax=min(r);
syms t
y=x1+t*Pk;
ft(t)=subs(f,x,y);
[r1]=find0618(ft,0,double(rmax),0.00001);
x1=x1+r1*Pk;
break;
end
end
count=count+1;
if isempty(M)
break;
end
if min(u)>=0
break;
end
end
X=x1;
FMIN=subs(f,x,X);
K=count;
end
function [A1,A2,b1,b2,k ] = fenjie( A,b,x )
% 分解起作用约束
A=A;
b=b;
x0=x;
k=0;q=0;
s=size(A);
A1=zeros(s(1),s(2));
A2=zeros(s(1),s(2));
b1=zeros(s(1),1);
b2=zeros(s(1),1);
for i=1:s(1)
gi=A(i,:)*x0-b(i);
if abs(gi)<0.0000001
k=k+1;
A1(k,:)=A(i,:);
b1(k,1)=b(i);
else
q=q+1;
A2(q,:)=A(i,:);
b2(q,1)=b(i);
end
end
if k>0
A1=A1(1:k,:);
b1=b1(1:k,:);
else
A1=[];
end
if q>0
A2=A2(1:q,:);
b2=b2(1:q,:);
end
end
(2) 数值算例:
Min f(x)= 2*x1^2+2*x2^2 – 2*x1*x2 – 4*x1 – 6*x2
S.T. –x1 – x2>= –2
–x1 – 5*x2>= –5
x1>=0
x2>=0
初值x0=[0;0]
在matlab command window里输入
syms x1 x2
f=2*x1^2+2*x2^2-2*x1*x2-4*x1-6*x2;
A=[-1 -1;-1 -5;1 0;0 1];
b=[-2;-5;0;0];
x1=[0;0];
x=[x1;x2];
e=0.01;
[X,FMIN,K]=tidutouying(f,A,b,x1,x,e)
X =
1.1292
0.7742
FMIN =
-7.1613
N =
16
