Zoutendijk法应用实例

问题：
在进行新装备设计开发时，必须考虑到装备可靠性方面的要求。

所谓可靠性，指的是装备在规定的时间和工作条件下，无故障的完成任务的能力。

装备的组成形式可以分为串联和并联两种的情况。

串联结构如图1所示，可靠性计算方式为；并联结构如图2所示，可靠性计算方式为，分别表示系统、部件1和部件2的可靠性。

为了提高系统的总体可靠性，我们希望单个部件的可靠性水平越高越好。

但部件的可靠性水平与成本密切关联，可靠性水平越高，成本也就越高。

因此，我们希望在一定的成本约束下，使系统的总体可靠性水平达到最高。

12
图1 串联结构
1
2
图2 并联结构
现有一装备系统，其组成结构如图3所示。

可以看出，该系统由5个部件组成。

假定各部件的基础可靠性为，对应的基础成本为，可靠性每提高0.01，部件成本提升量为，即可靠性增长量与成本增长量的比例系数。

具体值参见表1。

2
3
4
5
1
图3 装备结构图
表1 部件性能参数
部件编号基础可靠性基础成本成本提升量(可靠性每增加0.01)
1 0.90 10 1
2 0.92 18 1.5
3 0.95 20 2
4 0.89 9 0.9
5 0.95 20 2
现要求在总成本不超过90的前提，使得装备的总体可靠性水平达到最高。

提示：可靠性是在[0,1]区间内的连续量；总体成本为各部件成本之和。

分析：
设5个部件的可靠性提高量分别为
1
R
∆、
2
R
∆、
3
R
∆、
4
R
∆、
5
R
∆
则系统可靠性
5
4
3
2
1
))
1(
)
1(
)
1(
1(R
R
R
R
R
R
s
⋅
-
⋅
-
⋅
-
-
⋅
=
)95.0())11.0)(05.0)(08.0(1()90.0(54321R R R R R ∆+⋅-∆-∆-∆+⋅∆+= 约束条件为：
∈∆+190.0R [0 1]、 ∈∆+292.0R [0 1]、 ∈∆+395.0R [0 1]、 ∈∆+489.0R [0 1]、 ∈∆+595.0R [0 1]；
054321≥∆∆∆∆∆R R R R R 、、、、；
90)29.025.1(10020920181054321≤∆+∆+∆+∆+∆+++++R R R R R
将其化简后得
⎪⎪⎪
⎪⎪⎪
⎪⎪⎪⎩
⎪
⎪⎪
⎪⎪⎪
⎪⎪⎪⎨⎧≤∆-≤∆-≤∆-≤∆-≤∆-≤∆+∆+∆+∆+∆≤∆≤∆≤∆≤∆≤∆→∆+⋅-∆-∆-∆+⋅∆+-0
000013
.029.025.105
.011.005.008.01.0:
min )95.0())11.0)(05.0)(08.0(1()90.0(543215
4321
5432154321R R R R R R R R R R R R R R R st R R R R R 由于问题计算量过大，借助于matlab 工具实现 源代码：
A=[1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 1;1 1.5 2 0.9 2;-1 0 0 0 0;0 -1 0 0 0;0 0 -1 0 0; 0 0 0 -1 0;0 0 0 0 -1];
b=[0.1 0.08 0.05 0.11 0.05 0.13 0 0 0 0 0]; x0=[0 0 0 0 0]'; Aeq=[]; beq=[];
format long; func;
sz=length(x0); [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);
%求解线性规划问题得到可行下降方向d0
s=diff_val(x0);%自定义函数diff_val(x0)作用是求所给函数在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(diff_val(x0),A1,p,Aeq,q,lb,ub);
d0;mn;
df=abs(s*d0);
if df<1e-14
'最优解d为：';
x0;
f_val=fval(x0);
k;
return
else
%进行一维搜索，求f(xk(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<vmax
vmax=v;
end
end
end
end
end
vmax;
h=fmin(x0,d0,vmax);
a=x0+h*d0;
f_val=fval(a);
x0=x0+h*d0;
end
回调函数：
●func 记录函数及其自变量信息
syms x1 x2 x3 x4 x5;
f=-(0.9+x1)*(0.95+x5)-(0.9+x1)*(0.95+x5)*(x2-0.08)*(x3-0.05)*(x4-0.11); x=[x1,x2,x3,x4,x5];
●fval 计算函数在x0处得函数值
function f_val=fval(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;
计算结果： x0 =
0.068373258286964 0.000000005045114 0.000000006952478 0.068373209985112 0.000000008008447
f_val =
-0.919801424130391 故当
8008447
0.9500000099851120.95837320 6952478
0.95000000 5045114
0.92000000 82869640.9683732554321=====R R R R R 时，得到最优值：41303910.91980142=s R。