约束最优化问题

约束最优化问题

一实习目的

1.熟练掌握科学与工程计算中常用的基本算法；

2.掌握分析问题，设计算法的能力；

3.掌握模块化程序设计的基本思想，注重模块的“高内聚，低耦合”；

4.采用自顶向下，逐步细化的编程思想完成程序书写；

5.牢固建立“清晰第一，效率第二”的软件设计观念；

6.掌握软件调试，测试的基本技能和方法；

7.提高科技报告的书写质量；

8.在掌握无约束最优化问题求解方法的前提下，对一般情形下的约束最优化问题进行研究，通过实习掌握外点罚函数法、内点罚函数法、乘子法、线性近似规划法和序列二次规划法在求解一般情形下的约束最优化问题的应用。

二问题定义及题目分析

问题1：

要求用外点罚函数法和内点罚函数法解决约束问题：

Min f(x)=错误！未找到引用源。

s.t. 错误！未找到引用源。

错误！未找到引用源。

错误！未找到引用源。

问题2：

要求用乘子法解决约束问题：

Min 错误！未找到引用源。

s.t. 错误！未找到引用源。

错误！未找到引用源。

(错误！未找到引用源。)

问题3：

要求用线性近似规划法和序列二次规划法解决约束问题:

Min 错误！未找到引用源。

s.t. 错误！未找到引用源。

错误！未找到引用源。

错误！未找到引用源。

错误！未找到引用源。

三程序概要设计

1.外点罚函数法

Step1. 给定初始点错误！未找到引用源。,罚参数序列{错误！未找到引用源。}(常取错误！未找到引用源。),精度错误！未找到引用源。，并令k=0；

Step2. 构造增广目标函数错误！未找到引用源。;

Step3. 求解无约束优化问题min 错误！未找到引用源。,x错误！未找到引用源。,其解记为错误！未找到引用源。；

Step4. （终止准则：惩罚项充分小，或等价地错误！未找到引用源。近似可行）若错误！未找到引用源。，或者错误！未找到引用源。,错误！

未找到引用源。,则得解错误！未找到引用源。，否则令k=k+1，转 Step2.

2.内点罚函数法：

Step1. 给定初始可行解错误！未找到引用源。，罚参数序列{错误！未找到引用源。}(常取错误！未找到引用源。),精度错误！未找到引用源。，并令

k=0；

Step2. 构造增广目标函数错误！未找到引用源。;

Step3. 求解无约束优化问题min 错误！未找到引用源。,x错误！未找到引用源。,其解记为错误！未找到引用源。；

Step4. （终止准则）若错误！未找到引用源。，则得解错误！未找到引用源。，否则令k=k+1，转 Step2.

3.乘子法：

Step1. 给定初始点错误！未找到引用源。，初始lagrange乘子错误！未找到引用源。,i错误！未找到引用源。罚参数序列{错误！未找到引用源。},

精度错误！未找到引用源。，并令k=0；

Step2. 构造增广目标函数错误！未找到引用源。

Step3. 求解无约束优化问题min 错误！未找到引用源。，x错误！未找到引用源。,其解记为错误！未找到引用源。；

Step4. （终止准则）若错误！未找到引用源。,则得解错误！未找到引用源。，否则令

K=k+1,转Step2.

4.线性近似规划法：

Step1. 给定初始点错误！未找到引用源。，步长限制错误！未找到引用源。,缩小系数错误！未找到引用源。。精度错误！未找到引用源。,并令k=0；Step2. 求解线性规划问题：min 错误！未找到引用源。

S.t. 错误！未找到引用源。

错误！未找到引用源。

错误！未找到引用源。

其解记为错误！未找到引用源。.

Step3. 若错误！未找到引用源。是约束优化问题的可行解，则令错误！未找到引用源。,转Step4；否则，取错误！未找到引用源。j=1,...,n,转Step2；Step4. （终止准则）若错误！未找到引用源。,且满足错误！未找到引用源。,或者若错误！未找到引用源。,则得问题的近似解错误！未找到引用源。；

否则令错误！未找到引用源。,k=k+1,转Step2。

5.序列二次规划法：

Step1. 给定初始点错误！未找到引用源。,令k=0；

Step2. 若错误！未找到引用源。满足约束问题的k-T条件，停止计算，得到解错误！未找到引用源。；否则，转Step3；

Step3. 解二次规划问题 min 错误！未找到引用源。

s.t. 错误！未找到引用源。

错误！未找到引用源。

得解错误！未找到引用源。;

Step4. 令错误！未找到引用源。,转Step2.

四实验结果

第一题运行结果为：x=（1，0）

第二题运行结果为：x=（1，1）

第三题运行结果为：x=（2，2）

五结果分析及总结

经过理论运算，程序所得结果与理论值相同。通过这次实习，了解了非线性约束条件下最优化问题的几种求解方法：罚函数法，乘子法和线性近似规划法。了解了不同方法之间的优点和不足，如外点罚函数法不能保证运算过程中的点在可行域内，但是对初始条件要求不高等。

六源代码

1，2，3题通用函数：

function [x,minf] = minNT(f,x0,var,eps)

format long;

if nargin == 3

eps = 1.0e-6;

end

tol = 1;

x0 = transpose(x0);

gradf = jacobian(f,var);

jacf = jacobian(gradf,var);

while tol>eps

v = Funval(gradf,var,x0);

tol = norm(v);

pv = Funval(jacf,var,x0);

p = -inv(pv)*transpose(v);

p = double(p);

x1 = x0 + p;

x0 = x1;

end

x = x1;

minf = Funval(f,var,x);

format short;

function end

function fv = Funval(f,varvec,varval)

var = findsym(f);

varc = findsym(varvec);

s1 = length(var);

s2 = length(varc);

m =floor((s1-1)/3+1);

varv = zeros(1,m);

if s1 ~= s2

for i=0: ((s1-1)/3)

k = findstr(varc,var(3*i+1));

index = (k-1)/3;

varv(i+1) = varval(index+1);

end

fv = subs(f,var,varv);

else

fv = subs(f,varvec,varval);

end

function end

第一题：

function [x,minf] = minConPF(f,x0,g,h,c1,p,var,eps) format long;

if nargin == 7

eps = 1.0e-6;

end

k = 0;

FE = 0;

for i=1:length(h)

FE = FE + (h(i))^2;

end

x1 = transpose(x0);

x2 = inf;

while 1

M = c1*p;

FF = M*FE;

gx = Funval(g,var,x1);

gF = 0;

for i=1:length(g)

if gx(i)<0

gF = gF+M*(g(i)^2);

end

end

SumF = f + FF + gF;

[x2,minf] = minNT(SumF,transpose(x1),var);

if norm(x2 - x1)<=eps

x = x2;

break;

else

c1 = M;

x1 = x2;

end

end

minf = Funval(f,var,x);

format short;

第二题：

function [x,minf] = minFactor(f,x0,g,h,v,M,alpha,gama,var,eps) format long;

if nargin == 9

eps = 1.0e-4;

end

FE = 0;

for i=1:length(h)

FE = h(i)^2;

end

x1 = transpose(x0);

x2 = inf;

while 1

FF = M*FE;

Fh = v*h;

gF=0;

gx = Funval(g,var,x1);

for i=1:length(g)

if gx(i)>0

if gx(i)

gF=gF-v(i)*g(i)+M*(g(i)^2);

else

gF=gF-(v(i)^2)/M;

end

end

end

SumF = f + FF - Fh + gF;

[x2,minf] = minNT(SumF,transpose(x1),var);

Hx2 = Funval(h,var,x2);

Hx1 = Funval(h,var,x1);

if norm(Hx2) < eps

x = x2;

break;

else

if Hx2/Hx1 >= gama

M = alpha*M;

x1 = x2;

else

v = v - M*transpose(Hx2);

x1 = x2;

end

end

end

minf = Funval(f,var,x);

format short;

第三题：

function [x,minf] = minXXJS(f,g,X,alpha,sita,gama,beta,var,eps) if nargin == 8

eps =1.0e-4;

end

N=size(X);

n=N(2);

FX=zeros(1,n);

while 1

for i=1:n

FX(i)=Funval(f,var,X(:,i));

end

[XS,IX]=sort(FX);

Xsorted=X(:,IX);

px=sum(Xsorted(:,1:(n-1)),2)/(n-1);

Fpx=Funval(f,var,px);

SumF=0;

for i=1:n

SumF=SumF+(FX(IX(i))-Fpx)^2;

end

SumF=sqrt(SumF/n);

if SumF<=eps

x=Xsorted(:,1);

break;

else

bcon_1=1;

cof_alpha=alpha;

while bcon_1

x2=px+cof_alpha*(px-Xsorted(:,n)); gx2=Funval(g,var,x2);

if min(gx2)>=0

bcon_1=0;

else

cof_alpha=sqrt(cof_alpha);

end

end

fx2=Funval(f,var,x2);

if fx2

cof_gama=gama;

bcon_2=1;

while bcon_2

x3=px+cof_gama*(x2-px);

gx3=Funval(g,var,x3);

if min(gx3)>=0

bcon_2=0;

else

cof_gama=sqrt(cof_gama); end

end

fx3=Funval(f,var,x3);

if fx3

Xsorted(:,n)=x3;

X=Xsorted;

continue;

else

Xsorted(:,n)=x2;

X=Xsorted;

continue;

end

else

if fx2

Xsorted(:,n)=x2;

X=Xsorted;

continue;

else

if fx2

Xsorted(:,n)=x2;

end

cof_beta=beta;

bcon_3=1;

while bcon_3

x4=px+cof_beta*(Xsorted(:,n)-px);

gx4=Funval(g,var,x4);

if min(gx4)>=0

bcon_3=0;

else

cof_beta=cof_beta/2;

end

end

fx4=Funval(f,var,x4);

FNnew=Funval(f,var,Xsorted(:,n));

if fx4

Xsorted(:,n)=x4;

X=Xsorted;

continue;

else

x0=Xsorted(:,1);

for i=1:n

Xsorted(:,i)=x0+sita*(Xsorted(:,i)-x0); end

end

end

end

end

X=Xsorted;

end

minf=Funval(f,var,x);