DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx

1、DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码DEA算法学习笔记系列（三）一次性求解CCR模型所有DMU参数效率、规模效益、有效性特征、调整值的matlab代码1编写目的1.1Excel一次只能计算一个DMUDEA的CCR模型，他的对偶模型如下图： 很多人通过EXCEL提供的一个插件进行计算，如下图所示： 但是，这种方法有以下不足：（1）每次只能计算一个DMU，如果有多个DMU，那么需要人工重复计算过程多次；（2）通过Excel计算，只能得到，没法得到各个，所以，也无法直接判断是规模效益递增还是递减；（3）没发直接得到的值，也无法直接判断

2、DMU是弱DEA有效，还是DEA有效1.2Matlab编程一次性计算所有DMU的效率、有效性、调整值文章通过编写Matlab程序，实现一次性对所有DMU计算效率、有效性（根据以及所有的汇总值）、调整值（根据）。2Matlab求解线性规划2.1系统函数说明2.1.1 调用格式：x, fval, exitflag, output, lambda=linprog(f,A,b,Aeq,beq,lb,ub,x0,options) 注意：这个函数求的是最小值2.1.2输入参数说明输入f: 线性规划中的C，价值系数矩阵，一行n列的向量输入A,b:作有等式约束的问题。若没有不等式约束，则令A= 、b= 具体参

3、考例子输入lb ,ub：变量x的下界和上界输入x0：初值点输入options：为指定优化参数进行最小化Display 显示水平， 选择off 不显示输出；选择iter显示每一 步迭代过程的输出；选择final 显示最终结果。MaxFunEvals 函数评价的最大允许次数Maxiter 最大允许迭代次数TolX x处的终止容限2.1.3返回值说明返回值x：为最优解向量。返回值 fval：返回解x处的目标函数值exitflag ：描述函数计算的退出条件：若为正值，表示目标函数收敛于解x处；若为负值，表示目标函数不收敛；若为零值，表示已经达到函数评价或迭代的最大次数。1Function conver

4、ged to a solutionx.0Number of iterations exceededoptions.MaxIter.-2No feasible point was found.-3Problem is unbounded.-4NaNvalue was encountered during execution of the algorithm.-5Both primal and dual problems are infeasible.-7Search direction became too small. No further progress could be made.out

5、put 返回优化信息：output.iterations表示迭代次数；output.algorithm表示所采用的算法；outprt.funcCount表示函数评价次数。lambda 返回x处的拉格朗日乘子。它有以下属性： lambda.lower-lambda的下界； lambda.upper-lambda的上界； lambda.ineqlin-lambda的线性不等式； lambda.eqlin-lambda的线性等式。2.2简单例子，用代码求解2.2.1例1max f=-2x1+3x2s.t x1 + 2x2 8 4x1 16 4x2 12 x1,x2,x302.2.1.1求解过程（1）

6、先将目标函数转化成最小值问题：min(-f)= -2x1 - 3x2（2）计算代码： f=-2 -3; %相当于CA=1 2;4 0;0 4;b=8;16;12; x,fval=linprog(f,A,b) f=fval*(-1) 2.2.1.2结果x = 4.0000 2.0000fval = -14.0000f = 14.00002.2.2例2：minf=5x1-x2+2x3+3x4-8x5s.t2x1+x2-x3+x4-3x56 2x1+x2-x3+4x4+x57 0xj15j=1,2,3,4,5 2.2.2.1求解过程f=5 -1 2 3 -8; A=-2 1 -1 1 -3;2 1

7、-1 4 1; b=6;7; lb=0 0 0 0 0; ub=15 15 15 15 15; %注意：这个例题的xj除了大于0之外，还有小于15的条件x,fval=linprog(f,A,b,lb,ub) 2.2.2.2结果：x = 0.0000 0.0000 8.0000 0.000015.0000minf =-1042.2.3例题3（无解的例子） min f=5x1+x2+2x3+3x4+x5s.t2x1+x2-x3+x4-3x51 2x1+3x2-x3+2x4+x5-2 0xj1j=1,2,3,4,5程序： f=5 1 2 3 1; A=-2 1 -1 1 -3;2 3 -1 2 1;

8、 b=1;-2; lb=0 0 0 0 0; ub=1 1 1 1 1; x,fval,exitflag,output,lambda=linprog(f,A,b,lb,ub)运行结果： Exiting: One or more of the residuals, duality gap, or total relative error has grown 100000 times greater than its minimum value so far: the primal appears to be infeasible (and the dual unbounded). (The du

9、al residual TolFun=1.00e-008.)x = 0.0000 0.0000 1.1987 0.00000.0000fval =2.3975exitflag =-1output = iterations: 7 cgiterations: 0 algorithm: lipsollambda = ineqlin: 2x1 double eqlin: 0x1 double upper: 5x1 double lower: 5x1 double 显示的信息表明该问题无可行解。所给出的是对约束破坏最小的解。2.2.4例4（需要标准化的例子，一个等式的例子） 建立数学模型：max f=0

10、.15x1+0.1x2+0.08 x3+0.12 x4 s.tx1-x2- x3- x40 x2+ x3- x40 （需要标准化）x1+x2+x3+ x4=1 （一个等式的例子）xj0j=1,2,3,42.2.4.1求解过程（1）将其转换为标准形式：min z=-0.15x1-0.1x2-0.08 x3-0.12 x4 s.tx1-x2- x3- x40 -x2- x3+ x40 （需要标准化）x1+x2+x3+ x4=1xj0j=1,2,3,4（2）代码f = -0.15;-0.1;-0.08;-0.12; A =1 -1 -1 -1 0 -1 -1 1; b = 0; 0; Aeq=1 1

11、 1 1; %注意：等式在这里表示beq=1; lb = zeros(4,1); x,fval,exitflag = linprog(f,A,b,Aeq,beq,lb) f=-fval 2.2.4.2结果x = 0.5000 0.2500 0.0000 0.2500fval =-0.1300exitflag = 1f =0.1300 即4个项目的投资百分数分别为50%，25%，0,25%时可使该公司获得最大的收益，其最大收益可到达13%。过程正常收敛。2.2.5例5min f=-3x1+x2+x3s.t x1 - 2x2 +x3 11 -4x1 +x2 +2x3 3 -2x1 + x2 =1

12、x1,x2,x302.2.5.1求解过程（1）将其转换为标准形式：min f=-3x1+x2+x3s.t x1 - 2x2 +x3 11 4x1 -x2 -2x3 -3 -2x1 + x2 =1 x1,x2,x30（2）计算代码： f=-3 1 1; %相当于CA=1 -2 1;4 -1 -2;b=11;-3; Aeq=-2 0 1; %注意：等式在这里表示beq=1; x,fval=linprog(f,A,b,Aeq,beq) 2.2.5.2结果x = 4.0000 1.0000 9.0000fval = -2.00002.2.6例子6：松弛变量为基变量用等式重解例1max f=-2x1+3

13、x2s.t x1 + 2x2 8 4x1 16 4x2 12 x1,x2,x302.2.6.1求解过程（1）将其转换为标准形式，限制条件化为等式min f=-2x1-3x2+0x3+0x4+0x5s.t x1 + 2x2 +x3 =8 4x1 +x4 =16 4x2 +x5 =12 x1,x2,x30（2）计算代码： f=-2 -3 0 0 0; %相当于C，包含松弛变量A= ;b= ; Aeq=1 2 1 0 0;4 0 0 1 0;0 4 0 0 1;beq=8;16;12; lb = zeros(5,1); %5列1行的0，表示xi大等于0，不加这个条件，结果不一样x,fval=linp

14、rog(f,A,b,Aeq,beq,lb) f=fval*(-1) 2.2.6.2结果x = 4.0000 2.0000 0.0000 以下三个是松弛变量，例子1中没有计算出来 0.0000 4.0000fval = -14.0000f = 14.00002.3自定义Matlab函数求解线性规划（从excel读数据）2.3.1简单版：MyLinprog读取给定文件中数据，返回计算结果（1）数据格式如下：（2）函数代码function x,fval=MyLinprog(sDataPath)% 从Excel读取数据，计算线性规划，返回结果% 适用于计算最小值；都是等式（不等式要转换为等式）%-从e

15、xcel读取数据Aeq=xlsread(sDataPath,A)C=xlsread(sDataPath,C)beq=xlsread(sDataPath,b)lb=xlsread(sDataPath,x_LowLimit)x,fval=linprog(C,Aeq,beq,lb) （3）调用自定义函数，求解例子6的结果sDataPath=d:deaDEA_Data.xlsxa,b=MyLinprog(d:deaDEA_Data.xlsx)3DEA模型之CCR简介DEA第一个模型被命名为CCR模型。从生产函数的角度看,这一模型是用来研究具有多个输入,特别是具有多个输出的“生产部门”同时为“规模有效”

16、与“技术有效”的十分理想且卓有成效的方法。3.1CCR理论模型 例1： 某公司有甲、乙、丙三个企业，为评价这几个企业的生产效率，收集到反映其投入（固定资产年净值x1、流动资金x2、职工人数x3）和产出（总产值y1、利税总额y2）的有关数据如下表：企业指标甲乙丙x1（万元）41527x2（万元）1545x3（万元）825y1（万元）602224y2（万元）1268分析过程：对于第一个企业，产出综合值为60u1+12u2,投入综合值4v1+15v2+8v3，其中u1 u2 v1 v2 v3分别为产出与投入的权重系数。我们定义第一个企业的生产效率为：总产出与总投入的比：对其线性化过程如下：对应的对偶

17、问题如下（转换过程以及含义，参考文档DEA算法学系列之二：线性规划对偶问题与经济含义20171005V1.5）：4CCR模型计算过程一个决策单元的计算过程4.1例题说明例题8：以1997年全部独立核算企业为对象,对安徽、江西、湖南和湖北四省进行生产水平的比较。投入要素取固定资产净值年平均余额(亿元),流动资金年平均余额及从业人员(万人),产出要素取总产值(亿元)和利税总额(亿元).安徽江西湖南湖北固定资产932.66583.08936.841306.56流动资金980.45581.64849.311444.30从业人员401.8294.2443.20461.00利税总额179.2949.761

18、44.20181.41总产值2196.09930.221659.042662.21全要素相对生产率(即DEA评价值)1.0000.71400.92851.000排序13214.2基于理论构建模型湖南省4.3调整形式，以利于线性规划函数求解4.4按照自定义函数，构造excel文件4.4.1矩阵A的格式和说明表达的意思是：注意：E这列，对应的系数，因为这里分析的是湖南的，所以用的是湖南的数据；这个不需要事先找到一个可行基F列之后的是原来限制条件中的松弛变量的系数，注意最后两列4.4.2价值向量系数矩阵C的格式和说明总共11列，表示11个变量，顺序和A中的顺序要一致说明：e这列，对应的价值系数，表达

19、的意思就是：4.4.3资源限制矩阵b的格式和说明注意对比：4.4.4X取值条件的限制对应：说明：其中变量无限制，表示正负无穷大，这里用一个比较大的负数表示下限为无穷大4.5调用自定义函数（MyLinprog）求解指定决策单元模型sDataPath=d:deaDEA_Data.xlsxx,fval=MyLinprog(d:deaDEA_Data.xlsx)4.6计算结果评价4.6.1最优值说明：最优值的一般含义：在生产可能集T内，在产出量不变的情况下，尽量将投入量X0按同一比例减少.如果投入量X0不能按同一比例减少,即线性规划的最优值=1,则认为,决策单元j0既为技术有效,也为规模有效。表示：在

20、生产可能集T内，在产出量不变的情况下，投入量X能按同一比例减少，那么就可以认为：决策单元j0不为技术有效,或不为规模有效.4.6.2各变量的值意味着：4.6.3模型效率分析4.6.3.1判断规则：4.6.3.2湖南的判断结果因为：，不等于1，说明不是DEA有效；4.7调整方案在产出不变的情况，减少投入，或者在投入不变的情况下，增加产出5计算CCR模型的matlab函数所有决策单元在一个决策单元的基础上，根据CCR模型的特点，进行了改进，使能够一次性完成所有决策单元的计算。5.1程序代码（可直接运行）function DMUSita,lambda,S,AAdjusted,ADifferent=M

21、yDEA_ALL_DMU(sDataPath)%计算所有DMU的主函数%1、读取数据,数据要求，每一行表示一个指标，输入指标在上，输出指标在下；每一个列表示一个DMUIOData=xlsread(sDataPath,Input_Output);nParameter,nDMU=size(IOData);AData=IOData(:,2:nDMU);nDMU=nDMU-1; nInput=size(find(IOData(:,1)=1),1); DMUSita=;lambda=;S=AAdjusted=;%每一列代表一个决策单元DMU,每一行代表一个变量for iDMU=1:nDMU %计算第i个

22、DMU x,fval,exitflag=MyDEA_ONE_DMU(iDMU,AData,nInput,nParameter,nDMU); %记录效力值 lambda=lambda x(2:nParameter,:); S=S x(nDMU+2:(nDMU+2)+nParameter-1); tempSumlambda=sum(lambda(:,iDMU); tempSumSita=sum(S(:,iDMU); DMUSita=DMUSita fval;tempSumlambda;tempSumSita;exitflag;%3行n列，第一行是sita;%3行n列，第2行是lambda和；第三行

23、是松弛变量和 %调整 Input_Adjusted=AData(1:nInput,iDMU)*fval-x(nDMU+2:nInput+(nDMU+2)-1); Output_Adjusted=AData(nInput+1:nParameter,iDMU)+x(nInput+nDMU+2:nDMU+nParameter+1); DMUAdjusted=Input_Adjusted;Output_Adjusted; AAdjusted=AAdjusted DMUAdjusted;endADifferent=AAdjusted-AData; DMUSita(find(DMUSita0.000001

24、)=0;ADifferent(find(abs(ADifferent)0.000001)=0;AAdjusted(find(AAdjusted0.000001)=0;S(find(S0.000001)=0; function x,fval,exitflag=MyDEA_ONE_DMU(iDMU,AData,nInput,nParameter,nDMU)%准备数据%A,第一列是sita的系数；中间是lamnda的系数；最后是松弛变量的系数% 第一列，投入变量系数第i个决策单元投入的相反数，产出变量系数为0sitaA=-AData(1:nInput,iDMU);zeros(nParameter-n

25、Input,1);tempI=eye(nParameter)%产出变量的松弛变量系数为-1for i=nInput+1:nParameter tempI(:,i)=tempI(:,i)*-1;endAeq=sitaA AData tempI; C=1 zeros(1,nDMU) -0.00000001+zeros(1,nParameter);% 投入变量b为0；产出变量的b为第i个决策单元投入beq=zeros(nInput,1);AData(nInput+1:nParameter,iDMU); lb=-100000 zeros(1,nDMU+nParameter);ub=1 1000000+zeros(1,nDMU+nParameter);%ub=;x0=0 zeros(1,nDMU) beq;%options = optimset(Display,iter,MaxFunEvals ,1000,Maxiter,1000)options = optimset(LargeScale, off)x,fval,exitflag,output,lambda=linprog(C