求一DEA拓展模型（多投入单产出结构下的基于单要素投入最小化模型）的Matlab程序

原问：兄弟我实在忙不过来研究这个了，时间有限，matlab只会自己跑跑别人写好的程序，所以在此悬赏征集高人援手（倘若觉着工作量大，出价不够，价格可再议）！现已解决，现分享之~有所不对的地方请大家批评指正~
模型表述如下：

Energy use efficiency in U.S. manufacturing.pdf (243.55 KB)

2009-7-15 14:31:53 上传

具体模型见附文。

不论援手与否，看过此贴者，兄弟在此先谢过了~~

顶起~~ 研究中，望诸位能多少给点儿意见~~

像这样的操作，我建议楼主采用lingo编写，
过程简单易行一些
真的
1# xingjwchina

谢谢楼上的意见，问题已自行解决，经验证，结果是合适的。
现给出matlab程序代码，望大家提供意见，讨论、批评，欢迎之至~~

代码（1）： 多投入单产出结构下的基于单要素投入最小化模型


dt1=dt';
x=dt1(1:4,:);
y=dt1(5,:);
n=size(x',1);m=size(x,1);s=size(y,1);
eps=10^-10;
f=[zeros(1,n) -eps*ones(1,m+s) 1];
A=zeros(1,n+m+s+1);b=0;
LB=zeros(n+m+s+1,1);
LB(n+m+s+1)=-inf;UB=[];
c1=zeros(m+s,1);
aa=[0,0,0,0;0,0,0,0;0,0,1,0;0,0,0,0];
bb=[1,0,0,0;0,1,0,0;0,0,0,0;0,0,0,1];
ke=zeros(n,1)
for i=1:n;
    Aeq=[x eye(m) zeros(m,s) -aa*x(:,i)
        y zeros(s,m) -eye(s) zeros(s,1)];
    c=bb*x(:,i);
    beq=[c
        y(:,i)];
   w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
end;
w
lambada=w(1:n,:)
s_min=w(n+1:n+m,:)'
s_plus=w(n+m+1:n+m+s,:)'
Beta1=w(n+m+s+1,:)'

代码（2）： 多投入单产出结构下的基于单要素投入最小化的超效率模型

dt1=dt';
x=dt1(1:4,:);
y=dt1(5,:);
n=size(x',1);m=size(x,1);s=size(y,1);
eps=10^-10;
f=[zeros(1,n) -eps*ones(1,m+s) 1];
A=zeros(1,n+m+s+1);b=0;
LB=zeros(n+m+s+1,1);
LB(n+m+s+1)=-inf;UB=[];
c1=zeros(m+s,1);
aa=[0,0,0,0;0,0,0,0;0,0,1,0;0,0,0,0];
bb=[1,0,0,0;0,1,0,0;0,0,0,0;0,0,0,1];
ke=zeros(n,1)
for i=1:n;
    Aeq=[x eye(m) zeros(m,s) -aa*x(:,i)
        y zeros(s,m) -eye(s) zeros(s,1)];
    c=bb*x(:,i);
    beq=[c
        y(:,i)];
   w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
end;
ke=w(n+m+s+1,:)';
for i=1:n;
    i
    zz=w(n+m+s+1,i);
    if zz>0.9999999;
        if i>1&i<n;
        k=2;
        for j=1:10000000;
            xx=x(:,i);
            Aeq=[x(:,1:(i-1)) k*xx x(:,(i+1):n) eye(m) zeros(m,s) -k*aa*xx
            y zeros(s,m) -eye(s) zeros(s,1)];
        c=k*bb*xx;
        beq=[c
            y(:,i)];
        ww(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
        if ww(n+m+s+1,i)>0.9999999;
            ke(i,1)=k
            break;
        end;
        k=k-0.001;
        end;
        elseif i<2;
        k=2;
        for j=1:10000000;
            xx=x(:,i);
            Aeq=[k*xx x(:,2:n) eye(m) zeros(m,s) -k*aa*xx
            y zeros(s,m) -eye(s) zeros(s,1)];
        c=k*bb*xx;
        beq=[c
            y(:,i)];
        ww(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
        if ww(n+m+s+1,i)>0.9999999;
            ke(i,1)=k
            break;
        end;
        k=k-0.001;
        end;         
        else;
        k=2;
        for j=1:10000000;
            xx=x(:,i);
            Aeq=[x(:,1:n-1) k*xx eye(m) zeros(m,s) -k*aa*xx
            y zeros(s,m) -eye(s) zeros(s,1)];
        c=k*bb*xx;
        beq=[c
            y(:,i)];
        ww(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
        if ww(n+m+s+1,i)>0.9999999;
            ke(i,1)=k
            break;
        end;
        k=k-0.001;
        end;
    end;
    end;
end;
w
lambada=w(1:n,:)
s_min=w(n+1:n+m,:)'
s_plus=w(n+m+1:n+m+s,:)'
Beta1=w(n+m+s+1,:)'
Beta2=ww(n+m+s+1,:)'
kex=ke

注：1、代码（2）在（1）上，进行了适当拓展，具体涵义参见唐启义老师的DPS说明书（很遗憾，他开发的DEA模型种类太少，只包括几种最基本模型，不然就考虑买了）；
        2、超效率模型的方法比较粗鲁，无非是将（1）中处于前沿单元的点外推之后，再反复减小外推比例，直到DEA有效，因此，计算量比较大，还请诸位提供高见！（但该方法已通过EMS得到验证，仅由于精度原因，在第4位小数点上有所误差）
        3、模型为自己论文所用，所以没有拓展到一般（只包含4个投入，一个产出），需要再根据使用要求进行调整。

自此，感叹！其实DEA编程也没有那么难，在论坛上收获很多，得到很大帮助~~~感谢版友、感谢版主、感谢人大经济论坛、感谢cctv~~~~~~~~~~

xingjwchina 发表于 2009-7-15 21:52
谢谢楼上的意见，问题已自行解决，经验证，结果是合适的。
现给出matlab程序代码，望大家提供意见，讨论、批评，欢迎之至~~

代码（1）： 多投入单产出结构下的基于单要素投入最小化模型


dt1=dt';
x=dt1(1:4,:);
y=dt1(5,:);
n=size(x',1);m=size(x,1);s=size(y,1);
eps=10^-10;
f=[zeros(1,n) -eps*ones(1,m+s) 1];
A=zeros(1,n+m+s+1);b=0;
LB=zeros(n+m+s+1,1);
LB(n+m+s+1)=-inf;UB=[];
c1=zeros(m+s,1);
aa=[0,0,0,0;0,0,0,0;0,0,1,0;0,0,0,0];
bb=[1,0,0,0;0,1,0,0;0,0,0,0;0,0,0,1];
ke=zeros(n,1)
for i=1:n;
    Aeq=[x eye(m) zeros(m,s) -aa*x(:,i)
        y zeros(s,m) -eye(s) zeros(s,1)];
    c=bb*x(:,i);
    beq=[c
        y(:,i)];
   w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
end;
w
lambada=w(1:n,:)
s_min=w(n+1:n+m,:)'
s_plus=w(n+m+1:n+m+s,:)'
Beta1=w(n+m+s+1,:)'

代码（2）： 多投入单产出结构下的基于单要素投入最小化的超效率模型

dt1=dt';
x=dt1(1:4,:);
y=dt1(5,:);
n=size(x',1);m=size(x,1);s=size(y,1);
eps=10^-10;
f=[zeros(1,n) -eps*ones(1,m+s) 1];
A=zeros(1,n+m+s+1);b=0;
LB=zeros(n+m+s+1,1);
LB(n+m+s+1)=-inf;UB=[];
c1=zeros(m+s,1);
aa=[0,0,0,0;0,0,0,0;0,0,1,0;0,0,0,0];
bb=[1,0,0,0;0,1,0,0;0,0,0,0;0,0,0,1];
ke=zeros(n,1)
for i=1:n;
    Aeq=[x eye(m) zeros(m,s) -aa*x(:,i)
        y zeros(s,m) -eye(s) zeros(s,1)];
    c=bb*x(:,i);
    beq=[c
        y(:,i)];
   w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
end;
ke=w(n+m+s+1,:)';
for i=1:n;
    i
    zz=w(n+m+s+1,i);
    if zz>0.9999999;
        if i>1&i0.9999999;
            ke(i,1)=k
            break;
        end;
        k=k-0.001;
        end;
        elseif i0.9999999;
            ke(i,1)=k
            break;
        end;
        k=k-0.001;
        end;         
        else;
        k=2;
        for j=1:10000000;
            xx=x(:,i);
            Aeq=[x(:,1:n-1) k*xx eye(m) zeros(m,s) -k*aa*xx
            y zeros(s,m) -eye(s) zeros(s,1)];
        c=k*bb*xx;
        beq=[c
            y(:,i)];
        ww(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
        if ww(n+m+s+1,i)>0.9999999;
            ke(i,1)=k
            break;
        end;
        k=k-0.001;
        end;
    end;
    end;
end;
w
lambada=w(1:n,:)
s_min=w(n+1:n+m,:)'
s_plus=w(n+m+1:n+m+s,:)'
Beta1=w(n+m+s+1,:)'
Beta2=ww(n+m+s+1,:)'
kex=ke

注：1、代码（2）在（1）上，进行了适当拓展，具体涵义参见唐启义老师的DPS说明书（很遗憾，他开发的DEA模型种类太少，只包括几种最基本模型，不然就考虑买了）；
        2、超效率模型的方法比较粗鲁，无非是将（1）中处于前沿单元的点外推之后，再反复减小外推比例，直到DEA有效，因此，计算量比较大，还请诸位提供高见！（但该方法已通过EMS得到验证，仅由于精度原因，在第4位小数点上有所误差）
        3、模型为自己论文所用，所以没有拓展到一般（只包含4个投入，一个产出），需要再根据使用要求进行调整。

自此，感叹！其实DEA编程也没有那么难，在论坛上收获很多，得到很大帮助~~~感谢版友、感谢版主、感谢人大经济论坛、感谢cctv~~~~~~~~~~

很大公无私呀，俺来学习学习

要是版主能给点儿奖励就好了~

学习啊，学习啊！！！！！111

很强！！！！！！！！！！！！1

好厉害啊~~~~我用到的都是最简单的矩阵求解的幼稚园程序

这个貌似不错，因为我以前用LINGO做的，是多投入，多产出的情况