fun.asv

用简单遗传算法SGA实现的关于控制图的程序

function ECT=fun(x,par)
n=x(1);
h1=x(2);
h2=x(3);
k1=x(4);
k2=x(5);
%   par=[M a b lmd1 lmd2 dlt1 dlt2 g T0 T1 T2 L0 L1 L2 alpha_U PL FJTJ];
M=par(1);
a=par(2); b=par(3);
lmd1=par(4);lmd2=par(5);
dlt1=par(6);dlt2=par(7);g=par(8);
T0=par(9);T1=par(10);T2=par(11);
L0=par(12);L1=par(13);L2=par(14);

alpha1=2*PHI(-k2);
alpha2=2*PHI(-k1);

beta11=PHI(k2-sqrt(n)*dlt1);
beta12=PHI(k1-sqrt(n)*dlt1);
beta21=PHI(k2-sqrt(n)*dlt2);
beta22=PHI(k1-sqrt(n)*dlt2);

r1=1-exp(-lmd1*h1);
r2=1-exp(-lmd2*h2);
r=1+exp(-lmd1*h1)*lmd2/(lmd1-lmd2)-exp(-lmd2*h1)*lmd1/(lmd1-lmd2);
r3=1-exp(-lmd1*h2);
r4=1+exp(-lmd1*h2)*lmd2/(lmd1-lmd2)-exp(-lmd2*h2)*lmd1/(lmd1-lmd2);
r5=1-exp(-lmd2*h1);
r6=1-exp(-lmd2*(g*n*T1));

% P的第1列
p11=(1-r1)*(1-alpha1);
p21=(1-r3)*(1-alpha1);
p31=(1-r1)*(1-alpha1);
p41=0;
p51=0;
p61=(1-r1)*(1-alpha1);
p71=0;
p81=0;
p91=(1-r1)*(1-alpha1);
% P的第2列
p12=(1-r1)*(alpha1-alpha2);
p22=(1-r3)*(alpha1-alpha2);
p32=(1-r1)*(alpha1-alpha2);
p42=0;
p52=0;
p62=(1-r1)*(alpha1-alpha2);
p72=0;
p82=0;
p92=(1-r1)*(alpha1-alpha2);
% P的第3列
p13=(1-r1)*alpha2;
p23=(1-r3)*alpha2;
p33=(1-r1)*alpha2;
p43=0;
p53=0;
p63=(1-r1)*alpha2;
p73=0;
p83=0;
p93=(1-r1)*alpha2;
% P的第4列
p14=(r1-r)*beta11;
p24=(r3-r4)*beta11;
p34=(r1-r)*beta11;
p44=(1-r5)*beta11;
p54=(1-r2)*beta11;
p64=(r1-r)*beta11;
p74=0;
p84=0;
p94=(r1-r)*beta11;
% P的第5列
p15=(r1-r)*(beta12-beta11);
p25=(r3-r4)*(beta12-beta11);
p35=(r1-r)*(beta12-beta11);
p45=(1-r5)*(beta12-beta11);
p55=(1-r2)*(beta12-beta11);
p65=(r1-r)*(beta12-beta11);
p75=0;
p85=0;
p95=(r1-r)*(beta12-beta11);
% P的第6列
p16=(r1-r)*(1-beta12);
p26=(r3-r4)*(1-beta12);
p36=(r1-r)*(1-beta12);
p46=(1-r5)*(1-beta12);
p56=(1-r2)*(1-beta12);
p66=(r1-r)*(1-beta12);
p76=0;
p86=0;
p96=(r1-r)*(1-beta12);
% P的第7列
p17=r*beta21;
p27=r4*beta21;
p37=r*beta21;
p47=r5*beta21;
p57=r2*beta21;
p67=r*beta21;
p77=beta21;
p87=beta21;
p97=r*beta21;
% P的第8列
p18=r*(beta22-beta21);
p28=r4*(beta22-beta21);
p38=r*(beta22-beta21);
p48=r5*(beta22-beta21);
p58=r2*(beta22-beta21);
p68=r*(beta22-beta21);
p78=(beta22-beta21);
p88=(beta22-beta21);
p98=r*(beta22-beta21);
% P的第9列
p19=r*(1-beta22);
p29=r4*(1-beta22);
p39=r*(1-beta22);
p49=r5*(1-beta22);
p59=r2*(1-beta22);
p69=r*(1-beta22);
p79=(1-beta22);
p89=(1-beta22);
p99=r*(1-beta22);

p=[ p11 p12 p13 p14 p15 p16 p17 p18 p19
    p21 p22 p23 p24 p25 p26 p27 p28 p29
    p31 p32 p33 p34 p35 p36 p37 p38 p39
    p41 p42 p43 p44 p45 p46 p47 p48 p49
    p51 p52 p53 p54 p55 p56 p57 p58 p59
    p61 p62 p63 p64 p65 p66 p67 p68 p69
    p71 p72 p73 p74 p75 p76 p77 p78 p79
    p81 p82 p83 p84 p85 p86 p87 p88 p89
    p91 p92 p93 p94 p95 p96 p97 p98 p99];
if min(p(:))<0
    disp('出错了：概率值P有负值！'); 
    ppos=find(p<0);
    p(ppos)=0;
end
%
% syms pi00 pi01 pi02 pi10 pi11 pi12 pi20 pi21 pi22
eq1=strcat('pi00*(',num2str(p(1,1)),'-1)+ pi01*',num2str(p(2,1)),'+ pi02*',num2str(p(3,1)),'+ pi10*',num2str(p(4,1)),'+ pi11*',num2str(p(5,1)),'+ pi12*',num2str(p(6,1)),'+ pi20*',num2str(p(7,1)),'+ pi21*',num2str(p(8,1)),'+ pi22*',num2str(p(9,1)),'=0');
eq2=strcat('pi00*',num2str(p(1,2)),'+ pi01*(',num2str(p(2,2)),'-1)+ pi02*',num2str(p(3,2)),'+ pi10*',num2str(p(4,2)),'+ pi11*',num2str(p(5,2)),'+ pi12*',num2str(p(6,2)),'+ pi20*',num2str(p(7,2)),'+ pi21*',num2str(p(8,2)),'+ pi22*',num2str(p(9,2)),'=0');
eq3=strcat('pi00*',num2str(p(1,3)),'+ pi01*',num2str(p(2,3)),'+ pi02*(',num2str(p(3,3)),'-1)+ pi10*',num2str(p(4,3)),'+ pi11*',num2str(p(5,3)),'+ pi12*',num2str(p(6,3)),'+ pi20*',num2str(p(7,3)),'+ pi21*',num2str(p(8,3)),'+ pi22*',num2str(p(9,3)),'=0');
eq4=strcat('pi00*',num2str(p(1,4)),'+ pi01*',num2str(p(2,4)),'+ pi02*',num2str(p(3,4)),'+ pi10*(',num2str(p(4,4)),'-1)+ pi11*',num2str(p(5,4)),'+ pi12*',num2str(p(6,4)),'+ pi20*',num2str(p(7,4)),'+ pi21*',num2str(p(8,4)),'+ pi22*',num2str(p(9,4)),'=0');
eq5=strcat('pi00*',num2str(p(1,5)),'+ pi01*',num2str(p(2,5)),'+ pi02*',num2str(p(3,5)),'+ pi10*',num2str(p(4,5)),'+ pi11*(',num2str(p(5,5)),'-1)+ pi12*',num2str(p(6,5)),'+ pi20*',num2str(p(7,5)),'+ pi21*',num2str(p(8,5)),'+ pi22*',num2str(p(9,5)),'=0');
eq6=strcat('pi00*',num2str(p(1,6)),'+ pi01*',num2str(p(2,6)),'+ pi02*',num2str(p(3,6)),'+ pi10*',num2str(p(4,6)),'+ pi11*',num2str(p(5,6)),'+ pi12*(',num2str(p(6,6)),'-1)+ pi20*',num2str(p(7,6)),'+ pi21*',num2str(p(8,6)),'+ pi22*',num2str(p(9,6)),'=0');
eq7=strcat('pi00*',num2str(p(1,7)),'+ pi01*',num2str(p(2,7)),'+ pi02*',num2str(p(3,7)),'+ pi10*',num2str(p(4,7)),'+ pi11*',num2str(p(5,7)),'+ pi12*',num2str(p(6,7)),'+ pi20*(',num2str(p(7,7)),'-1)+ pi21*',num2str(p(8,7)),'+ pi22*',num2str(p(9,7)),'=0');
eq8=strcat('pi00*',num2str(p(1,8)),'+ pi01*',num2str(p(2,8)),'+ pi02*',num2str(p(3,8)),'+ pi10*',num2str(p(4,8)),'+ pi11*',num2str(p(5,8)),'+ pi12*',num2str(p(6,8)),'+ pi20*',num2str(p(7,8)),'+ pi21*(',num2str(p(8,8)),'-1)+ pi22*',num2str(p(9,8)),'=0');
eq9=strcat('pi00*',num2str(p(1,9)),'+ pi01*',num2str(p(2,9)),'+ pi02*',num2str(p(3,9)),'+ pi10*',num2str(p(4,9)),'+ pi11*',num2str(p(5,9)),'+ pi12*',num2str(p(6,9)),'+ pi20*',num2str(p(7,9)),'+ pi21*',num2str(p(8,9)),'+ pi22*(',num2str(p(9,9)),'-1)=0');
eq10='pi00+ pi01+ pi02+ pi10+ pi11+ pi12+ pi20+ pi21+ pi22=1';


% eq0=strcat('pi00*(',num2str(p00),'-1)+ pi01*',num2str(p10),'+ pi02*',num2str(p20),'+ pi10*',num2str(p30),'+ pi11*',num2str(p40),'+ pi12*',num2str(p50),'+ pi20*',num2str(p60),'+ pi21*',num2str(p70),'+ pi22*',num2str(p80),'=0');
% eq1=strcat('pi00*',num2str(p01),'+ pi01*(',num2str(p11),'-1)+ pi02*',num2str(p21),'+ pi10*',num2str(p31),'+ pi11*',num2str(p41),'+ pi12*',num2str(p51),'+ pi20*',num2str(p61),'+ pi21*',num2str(p71),'+ pi22*',num2str(p81),'=0');
% eq2=strcat('pi00*',num2str(p02),'+ pi01*',num2str(p12),'+ pi02*(',num2str(p22),'-1)+ pi10*',num2str(p32),'+ pi11*',num2str(p42),'+ pi12*',num2str(p52),'+ pi20*',num2str(p62),'+ pi21*',num2str(p72),'+ pi22*',num2str(p82),'=0');
% eq3=strcat('pi00*',num2str(p03),'+ pi01*',num2str(p13),'+ pi02*',num2str(p23),'+ pi10*(',num2str(p33),'-1)+ pi11*',num2str(p43),'+ pi12*',num2str(p53),'+ pi20*',num2str(p63),'+ pi21*',num2str(p73),'+ pi22*',num2str(p83),'=0');
% eq4=strcat('pi00*',num2str(p04),'+ pi01*',num2str(p14),'+ pi02*',num2str(p24),'+ pi10*',num2str(p34),'+ pi11*(',num2str(p44),'-1)+ pi12*',num2str(p54),'+ pi20*',num2str(p64),'+ pi21*',num2str(p74),'+ pi22*',num2str(p84),'=0');
% eq5=strcat('pi00*',num2str(p05),'+ pi01*',num2str(p15),'+ pi02*',num2str(p25),'+ pi10*',num2str(p35),'+ pi11*',num2str(p45),'+ pi12*(',num2str(p55),'-1)+ pi20*',num2str(p65),'+ pi21*',num2str(p75),'+ pi22*',num2str(p85),'=0');
% eq6=strcat('pi00*',num2str(p06),'+ pi01*',num2str(p16),'+ pi02*',num2str(p26),'+ pi10*',num2str(p36),'+ pi11*',num2str(p46),'+ pi12*',num2str(p56),'+ pi20*(',num2str(p66),'-1)+ pi21*',num2str(p76),'+ pi22*',num2str(p86),'=0');
% eq7=strcat('pi00*',num2str(p07),'+ pi01*',num2str(p17),'+ pi02*',num2str(p27),'+ pi10*',num2str(p37),'+ pi11*',num2str(p47),'+ pi12*',num2str(p57),'+ pi20*',num2str(p67),'+ pi21*(',num2str(p77),'-1)+ pi22*',num2str(p87),'=0');
% eq8=strcat('pi00*',num2str(p08),'+ pi01*',num2str(p18),'+ pi02*',num2str(p28),'+ pi10*',num2str(p38),'+ pi11*',num2str(p48),'+ pi12*',num2str(p58),'+ pi20*',num2str(p68),'+ pi21*',num2str(p78),'+ pi22*(',num2str(p88),'-1)=0');
% eq9='pi00+ pi01+ pi02+ pi10+ pi11+ pi12+ pi20+ pi21+ pi22=1';

[pi00,pi01,pi02,pi10,pi11,pi12,pi20,pi21,pi22]=solve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq10,'pi00,pi01,pi02,pi10,pi11,pi12,pi20,pi21,pi22');

M_pi=[pi00 pi01 pi02;
    pi10 pi11 pi12;
    pi20 pi21 pi22;];

    
M_pi=double(M_pi);
if min(M_pi(:))<0
    disp('Pi(ij)<0');
    M_pi=1/9.0*ones(3,3);
end
T00=h1;
T01=h2;
T02=g*n+T0+h1;
T10=h1;
T11=h2;
T12=g*n+T1+h1;
T20=h1;
T21=h2;
T22=g*n+T2+h1;

M_T=[T00 T01 T02;
    T10 T11 T12;
    T20 T21 T22;];

C00=a*n+b+M*(h1-(lmd1+lmd2)/(lmd1*lmd2)-exp(-lmd1*h1)*lmd2/(lmd1*(lmd1-lmd2))+exp(-lmd2*h1)*lmd1/(lmd2*(lmd1-lmd2)));
C01=a*n+b+M*(h2-(lmd1+lmd2)/(lmd1*lmd2)-exp(-lmd1*h2)*lmd2/(lmd1*(lmd1-lmd2))+exp(-lmd2*h2)*lmd1/(lmd2*(lmd1-lmd2)));
C02=a*n+b+L0+M*(g*n+T0-(lmd1+lmd2)/(lmd1*lmd2)-exp(-lmd1*(g*n*T0))*lmd2/(lmd1*(lmd1-lmd2))+exp(-lmd2*(g*n+T0))*lmd1/(lmd2*(lmd1-lmd2))) + M*(h1-(lmd1+lmd2)/(lmd1*lmd2)-exp(-lmd1*h1)*lmd2/(lmd1*(lmd1-lmd2))+exp(-lmd2*h1)*lmd1/(lmd2*(lmd1-lmd2)));
C10=a*n+b+M*(h1-r5/lmd2);
C11=a*n+b+M*(h2-r2/lmd2);
C12=a*n+b+L1+M*(g*n+T1-r6/lmd2)+ M*(h1-(lmd1+lmd2)/(lmd1*lmd2)-exp(-lmd1*h1)*lmd2/(lmd1*(lmd1-lmd2))+exp(-lmd2*h1)*lmd1/(lmd2*(lmd1-lmd2)));
C20=a*n+b+M*h1;
C21=a*n+b+M*h2;
C22=a*n+b+L2+M*(g*n+T2)+ M*(h1-(lmd1+lmd2)/(lmd1*lmd2)-exp(-lmd1*h1)*lmd2/(lmd1*(lmd1-lmd2))+exp(-lmd2*h1)*lmd1/(lmd2*(lmd1-lmd2)));

M_C=[C00 C01 C02;
    C10 C11 C12;
    C20 C21 C22;];

% 目标函数：
ECT=sum(sum(M_pi.*M_C))/sum(sum(M_pi.*M_T));