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

上传人:b****6 文档编号:5827506 上传时间:2023-01-01 格式:DOCX 页数:29 大小:328.34KB
下载 相关 举报
DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx_第1页
第1页 / 共29页
DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx_第2页
第2页 / 共29页
DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx_第3页
第3页 / 共29页
DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx_第4页
第4页 / 共29页
DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx_第5页
第5页 / 共29页
点击查看更多>>
下载资源
资源描述

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

《DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx》由会员分享,可在线阅读,更多相关《DEA算法学习系列之三一次性求解CCR模型所有DMU参数效率规模效益有效性特征调整值的matlab代码.docx(29页珍藏版)》请在冰豆网上搜索。

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

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

 

DEA算法学习笔记系列

 

(三)一次性求解CCR模型所有DMU参数——效率、规模效益、有效性特征、调整值的matlab代码

1编写目的

1.1Excel一次只能计算一个DMU

DEA的CCR模型,他的对偶模型如下图:

很多人通过EXCEL提供的一个插件进行计算,如下图所示:

但是,这种方法有以下不足:

(1)每次只能计算一个DMU,如果有多个DMU,那么需要人工重复计算过程多次;

(2)通过Excel计算,只能得到θ,没法得到各个

,所以,也无法直接判断是规模效益递增还是递减;

(3)没发直接得到

的值,也无法直接判断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=[]

具体参考例子

输入lb,ub:

变量x的下界和上界

输入x0:

初值点

输入options:

为指定优化参数进行最小化

Display  显示水平,选择’off’不显示输出;

选择’iter’显示每一步迭代过程的输出;

选择’final’显示最终结果。

MaxFunEvals函数评价的最大允许次数

Maxiter最大允许迭代次数

TolX  x处的终止容限      

2.1.3返回值说明

返回值x:

为最优解向量。

返回值fval:

返回解x处的目标函数值

exitflag:

描述函数计算的退出条件:

若为正值,表示目标函数收敛于解x处;

若为负值,表示目标函数不收敛;

若为零值,表示已经达到函数评价或迭代的最大次数。

1

Functionconvergedtoasolution x.

 

0

Numberofiterationsexceeded options.MaxIter.

 

-2

Nofeasiblepointwasfound.

 

-3

Problemisunbounded.

 

-4

NaN valuewasencounteredduringexecutionofthealgorithm.

 

-5

Bothprimalanddualproblemsareinfeasible.

 

-7

Searchdirectionbecametoosmall.Nofurtherprogresscouldbemade.

output返回优化信息:

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例1

              maxf=-2x1+3x2

              s.t x1+2x2≤8

      4x1≤16

       4x2≤12

       x1,x2,x3≥0

2.2.1.1求解过程

(1)先将目标函数转化成最小值问题:

min(-f)=-2x1-3x2

(2)计算代码:

f=[-2-3];%%相当于C

A=[12;40;04];

b=[8;16;12];

[x,fval]=linprog(f,A,b)

f=fval*(-1)

2.2.1.2结果

x=

4.0000

2.0000

 

fval=

-14.0000

 

f=

14.0000

2.2.2例2:

minf=5x1-x2+2x3+3x4-8x5

s.t  –2x1+x2-x3+x4-3x5≤6

    2x1+x2-x3+4x4+x5≤7

    0≤xj≤15  j=1,2,3,4,5

2.2.2.1求解过程

f=[5-123-8];

A=[-21-11-3;21-141];

b=[6;7];

lb=[00000];

ub=[1515151515];

%%注意:

这个例题的xj除了大于0之外,还有小于15的条件

[x,fval]=linprog(f,A,b,[],[],lb,ub)

2.2.2.2结果:

x=

          0.0000

          0.0000

          8.0000

          0.0000

            15.0000

minf=

  -104

2.2.3例题3(无解的例子)

          minf=5x1+x2+2x3+3x4+x5

s.t  –2x1+x2-x3+x4-3x5≤1

                2x1+3x2-x3+2x4+x5≤-2

                  0≤xj≤1  j=1,2,3,4,5

程序:

      f=[51231];

      A=[-21-11-3;23-121];

      b=[1;-2];

      lb=[00000];

      ub=[11111];

      [x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb,ub) 

运行结果:

        

      Exiting:

Oneormoreoftheresiduals,dualitygap,ortotalrelativeerror

        hasgrown100000timesgreaterthanitsminimumvaluesofar:

        theprimalappearstobeinfeasible(andthedualunbounded).

        (Thedualresidual

x=0.0000

                  0.0000

                  1.1987

                  0.0000

                    0.0000

fval=

                  2.3975

exitflag=

                  -1

output=

          iterations:

7

          cgiterations:

0

        algorithm:

'lipsol'

lambda=

                  ineqlin:

[2x1double]

                eqlin:

[0x1double]

                upper:

[5x1double]

                lower:

[5x1double]

      显示的信息表明该问题无可行解。

所给出的是对约束破坏最小的解。

2.2.4   例4(需要标准化的例子,一个等式的例子)

      建立数学模型:

maxf=0.15x1+0.1x2+0.08x3+0.12x4

      s.t  x1-x2-x3-x4≤0

        x2+x3-x4≥0(需要标准化)

              x1+x2+x3+x4=1(一个等式的例子)

              xj≥0  j=1,2,3,4    

2.2.4.1求解过程

(1)将其转换为标准形式:

minz=-0.15x1-0.1x2-0.08x3-0.12x4

      s.t  x1-x2-x3-x4≤0

        -x2-x3+x4≤0(需要标准化)

              x1+x2+x3+x4=1

              xj≥0  j=1,2,3,4

(2)代码

f=[-0.15;-0.1;-0.08;-0.12];

A=  [1-1-1-1

                   0-1-11];

b=[0;0];

Aeq=[1111];%%注意:

等式在这里表示

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.2500

fval=

                          -0.1300

exitflag=

                              1

f=

0.1300

      即4个项目的投资百分数分别为50%,25%,0,  25%时可使该公司获得最大的收益,其最大收益可到达13%。

过程正常收敛。

2.2.5例5

              minf=-3x1+x2+x3

              s.t x1-2x2+x3≤11

      -4x1+x2+2x3≥3

       -2x1+x2=1

       x1,x2,x3≥0

2.2.5.1求解过程

(1)将其转换为标准形式:

              minf=-3x1+x2+x3

              s.t x1-2x2+x3≤11

      4x1-x2-2x3≤-3

       -2x1+x2=1

       x1,x2,x3≥0

(2)计算代码:

f=[-311];%%相当于C

A=[1-21;4-1-2];

b=[11;-3];

Aeq=[-201];%%注意:

等式在这里表示

beq=[1];

[x,fval]=linprog(f,A,b,Aeq,beq)

2.2.5.2结果

x=

4.0000

1.0000

9.0000

 

fval=

-2.0000

2.2.6例子6:

松弛变量为基变量——用等式重解例1

              maxf=-2x1+3x2

              s.t x1+2x2≤8

      4x1≤16

       4x2≤12

       x1,x2,x3≥0

2.2.6.1求解过程

(1)将其转换为标准形式,限制条件化为等式

              minf=-2x1-3x2+0x3+0x4+0x5

              s.t x1+2x2+x3=8

      4x1+x4=16

       4x2+x5=12

       x1,x2,x3≥0

(2)计算代码:

f=[-2-3000];%%相当于C,包含松弛变量

A=[];

b=[];

Aeq=[12100;40010;04001];

beq=[8;16;12];

lb=zeros(5,1);%5列1行的0,表示xi大等于0,不加这个条件,结果不一样

[x,fval]=linprog(f,A,b,Aeq,beq,lb)

f=fval*(-1)

2.2.6.2结果

x=

4.0000

2.0000

0.0000以下三个是松弛变量,例子1中没有计算出来

0.0000

4.0000

 

fval=

-14.0000

 

f=

14.0000

2.3自定义Matlab函数求解线性规划(从excel读数据)

2.3.1简单版:

MyLinprog——读取给定文件中数据,返回计算结果

(1)数据格式如下:

(2)函数代码

function[x,fval]=MyLinprog(sDataPath)

%从Excel读取数据,计算线性规划,返回结果

%适用于计算最小值;都是等式(不等式要转换为等式)

%--从excel读取数据

[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:

\dea\DEA_Data.xlsx'

[a,b]=MyLinprog('d:

\dea\DEA_Data.xlsx')

3DEA模型之CCR简介

DEA第一个模型被命名为CCR模型。

从生产函数的角度看,这一模型是用来研究具有多个输入,特别是具有多个输出的“生产部门”同时为“规模有效”与“技术有效”的十分理想且卓有成效的方法。

3.1CCR理论模型

例1:

某公司有甲、乙、丙三个企业,为评价这几个企业的生产效率,收集到反映其投入(固定资产年净值x1、流动资金x2、职工人数x3)和产出(总产值y1、利税总额y2)的有关数据如下表:

   企业指标

x1(万元)

4

15

27

x2 (万元)

15

4

5

x3 (万元)

8

2

5

y1 (万元)

60

22

24

y2 (万元)

12

6

8

分析过程:

对于第一个企业,产出综合值为60u1+12u2,投入综合值4v1+15v2+8v3,其中u1u2v1v2v3分别为产出与投入的权重系数。

我们定义第一个企业的生产效率为:

总产出与总投入的比:

对其线性化过程如下:

对应的对偶问题如下(转换过程以及含义,参考文档《DEA算法学系列之二:

线性规划对偶问题与经济含义20171005V1.5》):

4CCR模型计算过程——一个决策单元的计算过程

4.1例题说明

例题8:

以1997年全部独立核算企业为对象,对安徽、江西、湖南和湖北四省进行生产水平的比较。

投入要素取固定资产净值年平均余额(亿元),流动资金年平均余额及从业人员(万人),产出要素取总产值(亿元)和利税总额(亿元).

安徽

江西

湖南

湖北

固定资产

932.66

583.08

936.84

1306.56

流动资金

980.45

581.64

849.31

1444.30

从业人员

401.8

294.2

443.20

461.00

利税总额

179.29

49.76

144.20

181.41

总产值

2196.09

930.22

1659.04

2662.21

全要素相对生产率

(即DEA评价值)

1.000

0.7140

0.9285

1.000

排序

1

3

2

1

4.2基于理论构建模型——湖南省

4.3调整形式,以利于线性规划函数求解

4.4按照自定义函数,构造excel文件

4.4.1矩阵A的格式和说明

表达的意思是:

注意:

E这列,对应

的系数,因为这里分析的是湖南的,所以用的是湖南的数据;

这个不需要事先找到一个可行基

F列之后的是原来限制条件中的松弛变量的系数,注意最后两列

4.4.2价值向量系数矩阵C的格式和说明

总共11列,表示11个变量,顺序和A中的顺序要一致

说明:

e这列,对应

的价值系数,表达的意思就是:

4.4.3资源限制矩阵b的格式和说明

注意对比:

4.4.4X取值条件的限制

对应:

说明:

其中变量无限制,表示正负无穷大,这里用一个比较大的负数表示下限为无穷大

4.5调用自定义函数(MyLinprog)求解指定决策单元模型

sDataPath='d:

\dea\DEA_Data.xlsx'

[x,fval]=MyLinprog('d:

\dea\DEA_Data.xlsx')

4.6计算结果评价

4.6.1最优值

说明:

最优值的一般含义:

在生产可能集T内,在产出量不变的情况下,尽量将投入量X0按同一比例θ减少.如果投入量X0不能按同一比例θ减少,即线性规划的最优值θ=1,则认为,决策单元j0既为技术有效,也为规模有效。

表示:

在生产可能集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]=MyDEA_ALL_DMU(sDataPath)

%计算所有DMU的主函数

%1、读取数据,数据要求,每一行表示一个指标,输入指标在上,输出指标在下;每一个列表示一个DMU

[IOData]=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,每一行代表一个变量

foriDMU=1:

nDMU

%计算第i个DMU

[x,fval,exitflag]=MyDEA_ONE_DMU(iDMU,AData,nInput,nParameter,nDMU);

%记录效力值

lambda=[lambdax(2:

nParameter,:

)];

S=[Sx(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和;第三行是松弛变量和

%调整

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=[AAdjustedDMUAdjusted];

end

ADifferent=AAdjusted-AData;

DMUSita(find(DMUSita<0.000001))=0;

ADifferent(find(abs(ADifferent)<0.000001))=0;

AAdjusted(find(AAdjusted<0.000001))=0;

S(find(S<0.000001))=0;

function[x,fval,exitflag]=MyDEA_ONE_DMU(iDMU,AData,nInput,nParameter,nDMU)

%准备数据

%A,第一列是sita的系数;中间是lamnda的系数;最后是松弛变量的系数

%第一列,投入变量系数第i个决策单元投入的相反数,产出变量系数为0

sitaA=[-AData(1:

nInput,iDMU);zeros(nParameter-nInput,1)];

tempI=eye(nParameter)

%产出变量的松弛变量系数为-1

fori=nInput+1:

nParameter

tempI(:

i)=tempI(:

i)*-1;

end

[Aeq]=[sitaAADatatempI];

[C]=[[1]zeros(1,nDMU)-0.00000001+zeros(1,nParameter)];

%投入变量b为0;产出变量的b为第i个决策单元投入

beq=[zeros(nInput,1);AData(nInput+1:

nParameter,iDMU)];

lb=[-100000zeros(1,nDMU+nParameter)];

ub=[11000000+zeros(1,nDMU+nParameter)];

%ub=[];

x0=[0zeros(1,nDMU)beq'];

%options=optimset('Display','iter','MaxFunEvals',1000,'Maxiter',1000)

options=optimset('LargeScale','off')

[x,fval,exitflag,output,lambda]=linprog(C

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置:首页 > PPT模板 > 中国风

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1