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