鏈式網(wǎng)絡(luò)DEA與matlab應(yīng)用-二階段附加投入DEA

鏈式網(wǎng)絡(luò)DEA與matlab應(yīng)用：二階段附加投入的DEA鏈式網(wǎng)絡(luò)DEA包括2階段、3階段等多階段，我們考慮2階段附加投入的DEA模型。網(wǎng)絡(luò)DEA模型的形式是我們把決策單元分解為K個子過程，每個子過程都有相應(yīng)的投入變量和輸出變量，同時還有中間變量，中間變量既是上一子系統(tǒng)輸出的變量，又是下一子系統(tǒng)輸入的變量，我們把各個過程依照流程的次序排列好，上一子過程的中間產(chǎn)出變量要由下一個過程進行消費，通過相關(guān)中間變量，緊密聯(lián)系決策單元各子過程。對應(yīng)的一般模型為，DEA模型如下：第一階段CCR模型：不僅考慮到要求保證第一階段的效率值小于1，同時考慮附加約束即第二階段的效率值閾值為1，計算第一階段的最佳效率，記為。所以，第一階段的效率值的約束條件為。根據(jù)Charnes-Cooper變換，將第一階段模型轉(zhuǎn)換為線性方程，整體效率可以看成是函數(shù)，同理根據(jù)Charnes-Cooper變換，clear;clc;M=xlsread('C:\MATLAB74\work\new','Sheet1','A2:G28');%第一階段3個投入M的1-2列，產(chǎn)出為M的3列；%第二階段的新投入為M的4-5列，產(chǎn)出為M的6-7列X1=M(:,1:2);%第一階段初始投入矩陣X1，列數(shù)為DMU的個數(shù)n，行數(shù)為投入種類mZ=M(:,3:3);%第一階段產(chǎn)出矩陣Z，列數(shù)為DMU的個數(shù)n，行數(shù)為產(chǎn)出種類DX2=M(:,4:5);%第二階段追加投入x2，列數(shù)為DMU的個數(shù)n，行數(shù)為附加投入種類HY=M(:,6:7);%第二階段最終產(chǎn)出矩陣Y，列數(shù)為DMU的個數(shù)n，行數(shù)為最終產(chǎn)出種類sn=size(X1,1);m=size(X1,2);D=size(Z,2);H=size(X2,2);s=size(Y,2);theta1=zeros(n,1);theta2=zeros(n,1);w=cell(n,1);v=w;Q=w;u=w;options=optimset('display','off');fori=1:n%%第一次線性規(guī)劃f=-[Z(i,:),zeros(1,(m+H+s))];A=[Z,-X1,zeros(n,(H+s));-Z,zeros(n,m),-X2,Y];b=zeros(2*n,1);Aeq=[zeros(1,D),X1(i,:),zeros(1,(H+s))];beq=1;lb=zeros(1,(D+m+H+s));ub=inf*ones(1,(D+m+H+s));[w1,feval]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);theta1(i)=-feval;%%第二次線性規(guī)劃theta0=linspace(0,theta1(i));THETA=zeros(1,100);forj=1:100f=-theta0(j)*[zeros(1,m),zeros(1,D),zeros(1,H),Y(i,:)];A=[-X1,Z,zeros(n,H),zeros(n,s);zeros(n,m),-Z,-X2,Y];b=zeros(2*n,1);Aeq=[zeros(1,m),Z(i,:),X2(i,:),zeros(1,s);-theta0(j)*X1(i,:),Z(i,:),zeros(1,H),zeros(1,s)];beq=[1;0];lb=zeros(1,(m+D+H+s));ub=inf*ones(1,(m+D+H+s));[w2,feval]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);THETA(j)=-feval;end%[~,j]=max(THETA);f=-theta0(j)*[zeros(1,m),zeros(1,D),zeros(1,H),Y(i,:)];A=[-X1,Z,zeros(n,H),zeros(n,s);zeros(n,m),-Z,-X2,Y];b=zeros(2*n,1);Aeq=[zeros(1,m),Z(i,:),X2(i,:),zeros(1,s);-theta0(j)*X1(i,:),Z(i,:),zeros(1,H),zeros(1,s)];beq=[1;0];lb=zeros(1,(m+D+H+s));ub=inf*ones(1,(m+D+H+s));[x,feval]=linprog(f,A,b,Aeq,beq,lb,ub,[],options);theta2(i)=-feval;v{i}=x(1:m);