zoutendijk 可行方向法的matlab实现

(一) 、基本思想是：

给定一个可行点x(k)之后，用某种方法确定一个改进的可行方向dk，然后沿方向dk，求解一个有约束的线搜索问题，得极小点?k，令x(k?1)?x(k)??kdk，如果x(k?1)不是最优解，则重复上述步骤。可行方向法就是利用线性规划方法来确定dk的。

1) 、线性约束问题：

设x是问题

?min f(x)? Ax?b,Ex?e ?s.t.? x?Rn?的一个可行解，假定A1x?b1，A2x?b2，其中

?A??b?A??1?，b??1?

?A2??b2?则一个非零向量d是在点x点的一个可行方向，当且仅当A1d?0，Ed?0；如果?Tf(x)d?0，则d是一改进方向。

2) 、非线性约束问题

设x是问题

?min f(x)? gi(x)?0, i?1,2,?,m ?s.t.?n x?R?的一个可行解，令S?x|x?Rn,gi(x)?0，I??x|gi(x)?0?，即I是x?S点紧约束的指标集，设f和gi(i?I)在x点可微，gi(i?I)在x点连续，如果

???Tf(x)d?0，?Tgi(x)d?0(i?I)，则d是一改进的可行方向。

(二) 、算法

1) 、线性不等式约束的Zoutendijk方法的计算步骤：

1.求一初始可行解x0。，令k＝1，转2。

TT 2.对于可行点xk，设A1xk?b1，A2xk?b2，AT?(A1T,A2)， bT?(b1T,b2)求解

?min z??f?x?Td?T?s.t A1d?0?fxdk=0，d问题P ，得最优解，如果计算结束， ??k1?? Ed?0? ?1?d?1,j?1,?,nj?xk是K—T点；否则转3。

3.求解线搜索问题

?min f(xk??dk) (a) ??s.t 0????max其中

?max??min ?b_d_|d_?0?,d_??0?? b_?b2?A2xk,d_?A2dk ????,d?0设?k为(a)式最优解，令xk?1?xk??kdk,k?k?1，返回2。

2) 、非线性不等式约束的Zoutendijk方法的计算步骤：

1) 选取允许误差?1?0，?2?0，求一初始可行点x(1)，令k?1，转2)。 2）确定指标集I(x(k))?i|gi(x(k))?0。

3) 若I(x(k))??，且?f(x(k))??1，计算结束，取x*?x(k)；若I(x(k))??，且?f(x(k))??1，令dk???f(x(k))，转6)；若I(x(k))??，转4)。

4) 令x?x(k)，求解线性规划问题(4-2)的最优解(zk,dk)；

??5) 若zk??2，计算结束，取x?x(k)；否则令dk＝dk，转6)。 6) 求出线搜索问题

?min f(x(k)??dk) ?s.t. 0????max?的最优解?k，其中?max?max?|x(k)??dk?S；令x(k?1)?x(k)??dk，

k?k?1，返回2)。

??(三) 、程序源码

1) 、主程序

简单说明：此程序可以处理线性和非线性问题，程序主要由label得值来判断，当label=1时运行线性约束部分，label=0时运行非线性约束部分 function [X0,f_val]=zoutendijk(A,b,x0,Aeq,beq,label)

%自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数 %自定义函数fval(x0)作用是求所给函数在x0出的函数值 format long; eps=1.0e-6;

x0=transpose(x0);%刚开始给的x0为行向量 func

sz=length(x0); if label==1 [m,n]=size(A);

%把A分解为A1,A2，其中A1为起作用约束 for k=1:1:100 A1=A; A2=A; b1=b; b2=b; for i=m:-1:1

if abs(A2(i,:)*x0-b2(i,:)) < 0.1 A2(i,:)=[];

b2(i,:)=[]; end end for i=m:-1:1

if abs(A1(i,:)*x0-b1(i,:))>=0.1 A1(i,:)=[]; b1(i,:)=[]; end end A1; A2; b1; b2;

i2=rank(A2); AE=[A1;Aeq]; [i1,j1]=size(AE); r=rank(AE); if r

'行不满秩' return end

if i2==0 '无效' return end

%求解线性规划问题得到可行下降方向d0 s=diff_val(x0); c=double(s); lb=-1*ones(sz,1); ub=ones(sz,1); k1=length(b1); k2=length(beq);

p=zeros(k1,1); q=zeros(k2,1);

[d0,mn,m1,m2,m3]=linprog(c,A1,p,Aeq,q,lb,ub); d0;mn;

df=abs(s*d0); if df<0.1 '最优解为：'

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' x0

f_val=fval(x0) k

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' return else

%进行一维搜索，求f(x(k+1))的最小值 b_=b2-A2*x0; d_=A2*d0; [dh,dl]=size(d_); ul=1; for i=1:1:dh if d_(i,:)>=0 u=1; else u=0; end ul=ul*u; end ul;b_;d_; vmax=inf;

if ul==0 vmax=inf; else

for i=1:1:dh if d_(i,:)>0

v=b_(i,:)/d_(i,:); if v

h=fmin(x0,d0,vmax); a=x0+h*d0; f_val=fval(a); x0=x0+h*d0; '****************' X0=x0 f_val=fval(x0) k

'****************' end end if label==0 for k=1:1:100 %确定指标集 'f(x)梯度：' sf=diff_val(x0) sf=eval(sf)

'f(x)梯度长度：' norm_s=norm(sf)

GL=length(G); G_copy=G; for i=1:1:GL

g=subs(G(i,:),x,x0); G(i,:)=g; end

G_zero=eval(G); for i=GL:-1:1

if abs(G_zero(i,:))>0.1 G_zero(i,:)=[]; G_copy(i,:)=[]; I=length(G_zero); end end

'x0时为零的g(x)：' G_copy '指标集I(x)：' I

add=-ones(I,1);

%根据指标集确定不同情况下的搜索方向 if I==0

if norm_s<=10 '最优解为：'

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' x0

f_val=fval(x0) k

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' return else

d0=-sf; %线搜索问题 vmax=100; d0=transpose(d0); h=fmin(x0,d0,vmax); x0=x0+h*d0; end else

%线性规划问题 grad=jacobian(G_copy,x); G_zero=subs(grad,x,x0); G_zero=[G_zero,add]; sf=[sf,-1];

'线性规划问题A矩阵：' Ac=[sf;G_zero] lb=-1*ones(sz,1); ub=ones(sz,1); p=zeros(I+1,1); c=zeros(1,sz); c=[c,1];

[dz,mn,m1,m2,m3]=linprog(c,Ac,p,[],[],lb,ub); dz; mn;

sd=length(dz); d1=dz(1:sd-1,1:1); z0=dz(sd,1); z0=abs(z0) if z0<0.01 '最优解为：'

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' x0

f_val=fval(x0)

k

'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' return else d0=d1; %线搜索问题 vmax=10000; h=fmin(x0,d0,vmax); a=x0+h*d0; x0=x0+h*d0; end end k end end

2) 、回调函数

? func 记录函数及其自变量信息，如： syms x1 x2;

f=2*x1^2+2*x2^2-2*x1*x2-4*x1-6*x2; x=[x1,x2];

? fval 计算函数在x0处得函数值 function f_val=fval(x0) x0=transpose(x0); func;

f_val=subs(f,x,x0);

? diff_val(x0)计算函数在x0处得导数值 function s=diff_val(x0) func

grad=jacobian(f,x); s=subs(grad,x,x0);

? fmin(x0,d0,vmax)求函数在[0,vmax]上的最小值 function h=fmin(x0,d0,vmax)

func syms h; a=x0+h*d0;

f_val=inline(subs(f,x,a)); if vmax==inf

min_h=fminbnd(f_val,0,10000); else

min_h=fminbnd(f_val,0,vmax); end

h=min_h;

(四) 、例题

? 线性问题（以例1为例说明数据输入及其最后结果）

例1

2?min f?x??2x12?2x2?2x1x2?4x1?6x2??s.t x1?x2?2? ? x1?5x2?5? ?x?01??? ?x2?0?3524?最优解x???,?,f?x????7.16

?3131?首先把func函数中相应例题的注释去掉； 然后在matlab中输入如下数据： clear clc

A=[1 1;1 5;-1 0;0 -1] b=[2 5 0 0]' x0=[0 0] Aeq=[] beq=[] label=1

最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label) 得到结果：

例2

2?min f?x??x12?4x2??s.t x1?x2?1?? 15x1?10x2?12 ? x?01??? x2?04?41?最优解x???,?,f?x???

5?55?A=[-1 -1;-15 -10;-1 0;0 -1]

b=[-1 -12 0 0]' x0=[0 2] Aeq=[] beq=[] label=1;

? 非线性问题 例

2?min f?x??x12?x1x2?2x2?6x1?2x2?12x3?2?15?s.t 2x12?x2 ?? -x1?2x2?x3?3? x,x,x?0123?

首先把func函数中相应例题的注释去掉； 然后在matlab中输入如下参数： clear clc A=[] b=[]

x0=[0 0 0] Aeq=[] beq=[] label=0

最后运行程序：zoutendijk(A,b,x0,Aeq,beq,label) 得到结果：

本文来源：https://www.bwwdw.com/article/kolf.html