381777078
2012-05-14, 16:05
function slove_Gini_Coefficient_2
%
%
%
%%
clc;
clear;
%%
n = 15;
m = 6;
W = 0.018 * 10^4;
Aeq = ones(1,n);
beq = W;
x0 = [0.0191 0.0191 0.0004 0.0089 0.0030 0.0057 0.0010 0.0018 0.0016 0.0010 0.0017 0.0046 0.0007 0.0012 0.0006]' * 10^4;
z = [10.070 47.602 6.974 158.565 532.687 1.557;
8.030 19.117 6.416 57.970 172.880 0.387 ;
1.711 2.659 0.274 3.069 58.347 0.036 ;
3.818 5.317 0.278 5.115 55.106 0.096 ;
4.410 7.216 0.370 14.322 44.301 0.120 ;
2.699 3.292 0.332 6.138 44.301 0.045 ;
3.818 3.038 0.228 3.069 35.657 0.069 ;
3.554 3.671 0.189 2.728 31.335 0.081 ;
4.673 6.457 0.145 10.571 28.093 0.084 ;
2.106 10.001 1.155 64.449 25.932 0.306 ;
6.977 4.431 0.172 6.479 16.208 0.081 ;
4.937 2.659 0.147 3.069 14.047 0.063 ;
2.304 3.671 0.135 2.387 10.805 0.027 ;
3.686 5.317 0.136 2.387 8.644 0.033 ;
3.028 2.152 0.089 0.682 2.161 0.015];
z(:,1) = z(:,1) * 10^4; z(:,2) = z(:,2) * 10^8;
z(:,3) = z(:,3) * 10^8; z(:,4) = z(:,4) * 10^8;
z(:,5) = z(:,5) * 10^4; z(:,6) = z(:,6) * 10^8;
P1 = 0.25;
P2 = 0.35; % 0.25
lb = x0 * (1 - P2).^4;
ub = x0 * (1 - P1).^4;
G0 = [0.4339 0.3781 0.3660 0.5003 0.3705 0.4445]';
options = optimset('Algorithm','interior-point','TolCon',1e-3);
[x f ] = fmincon(@(x)myfun(x,z),x0,[],[],Aeq,beq,lb,ub,@(x)mycon(x,z),options)
[w G] = get_w_G(x,z);
w
G
%%
function f = myfun(x,z)
% 目标函数
% z -> n *m (n = 15,m = 6)
%
[n m] = size(z);
y = zeros(n,m);
for i = 1 : n
y(i,:) = x(i) ./ z(i,:);
end
p = zeros(n,m);
sum_yj = sum(y);
for j = 1 : m
p(:,j) = y(:,j) / sum_yj(j);
end
e = zeros(m,1);
temp = 0.0;
for j = 1 : m
temp = 0.0;
for i = 1 : n
temp = temp + p(i,j) * log(p(i,j));
end
temp = - 1/(log(n)) * temp;
e(j) = temp;
end
w = zeros(m,1);
sum_ej = sum(1 - e);
for j = 1 : m
w(j) = (1 - e(j)) / sum_ej;
end
Z = [65.82 * 10^4 126.6 * 10^8 17.04 * 10^8 341 * 10^8 1.08 * 10^8 3 * 10^8]';
G = zeros(m,1);
temp = 0;
for j = 1 : m
for i = 1 : n
temp = temp + z(i,j) * (2 * sum(x(1:i)) - x(i)) / (Z(j) * 1.06 * 10^4);
end
G(j) = 1 - temp;
temp = 0.0;
end
temp = 0.0;
for j = 1 : m
temp = temp + w(j) * G(j);
end
f = temp;
function [c ceq] = mycon(x,z)
%
% 非线性约束条件
%
Z = [65.82 * 10^4 126.6 * 10^8 17.04 * 10^8 341 * 10^8 1.08 * 10^8 3 * 10^8]';
G0 = [0.4339 0.3781 0.3660 0.5003 0.3705 0.4445]';
m = 6;
n = 15;
G = zeros(m,1);
temp = 0;
for j = 1 : m
for i = 1 : n
temp = temp + z(i,j) * (2 * sum(x(1:i)) - x(i)) / (Z(j) * 0.018 * 10^4);
end
G(j) = 1 - temp;
temp = 0.0;
end
c = zeros(m,1);
for j = 1 : m
c(j) = G(j) - G0(j);
end
ceq = [];
%%
function [w G] = get_w_G(x,z)
%
%
[n m] = size(z);
y = zeros(n,m);
for i = 1 : n
y(i,:) = x(i) ./ z(i,:);
end
p = zeros(n,m);
sum_yj = sum(y);
for j = 1 : m
p(:,j) = y(:,j) / sum_yj(j);
end
e = zeros(m,1);
temp = 0.0;
for j = 1 : m
temp = 0.0;
for i = 1 : n
temp = temp + p(i,j) * log(p(i,j));
end
temp = - 1/(log(n)) * temp;
e(j) = temp;
end
w = zeros(m,1);
sum_ej = sum(1 - e);
for j = 1 : m
w(j) = (1 - e(j)) / sum_ej;
end
Z = [65.82 * 10^4 126.6 * 10^8 17.04 * 10^8 341 * 10^8 1.08 * 10^8 3 * 10^8]';
G = zeros(m,1);
temp = 0;
for j = 1 : m
for i = 1 : n
temp = temp + z(i,j) * (2 * sum(x(1:i)) - x(i)) / (Z(j) * 0.018 * 10^4);
end
G(j) = 1 - temp;
temp = 0.0;
end
%
%
%
%%
clc;
clear;
%%
n = 15;
m = 6;
W = 0.018 * 10^4;
Aeq = ones(1,n);
beq = W;
x0 = [0.0191 0.0191 0.0004 0.0089 0.0030 0.0057 0.0010 0.0018 0.0016 0.0010 0.0017 0.0046 0.0007 0.0012 0.0006]' * 10^4;
z = [10.070 47.602 6.974 158.565 532.687 1.557;
8.030 19.117 6.416 57.970 172.880 0.387 ;
1.711 2.659 0.274 3.069 58.347 0.036 ;
3.818 5.317 0.278 5.115 55.106 0.096 ;
4.410 7.216 0.370 14.322 44.301 0.120 ;
2.699 3.292 0.332 6.138 44.301 0.045 ;
3.818 3.038 0.228 3.069 35.657 0.069 ;
3.554 3.671 0.189 2.728 31.335 0.081 ;
4.673 6.457 0.145 10.571 28.093 0.084 ;
2.106 10.001 1.155 64.449 25.932 0.306 ;
6.977 4.431 0.172 6.479 16.208 0.081 ;
4.937 2.659 0.147 3.069 14.047 0.063 ;
2.304 3.671 0.135 2.387 10.805 0.027 ;
3.686 5.317 0.136 2.387 8.644 0.033 ;
3.028 2.152 0.089 0.682 2.161 0.015];
z(:,1) = z(:,1) * 10^4; z(:,2) = z(:,2) * 10^8;
z(:,3) = z(:,3) * 10^8; z(:,4) = z(:,4) * 10^8;
z(:,5) = z(:,5) * 10^4; z(:,6) = z(:,6) * 10^8;
P1 = 0.25;
P2 = 0.35; % 0.25
lb = x0 * (1 - P2).^4;
ub = x0 * (1 - P1).^4;
G0 = [0.4339 0.3781 0.3660 0.5003 0.3705 0.4445]';
options = optimset('Algorithm','interior-point','TolCon',1e-3);
[x f ] = fmincon(@(x)myfun(x,z),x0,[],[],Aeq,beq,lb,ub,@(x)mycon(x,z),options)
[w G] = get_w_G(x,z);
w
G
%%
function f = myfun(x,z)
% 目标函数
% z -> n *m (n = 15,m = 6)
%
[n m] = size(z);
y = zeros(n,m);
for i = 1 : n
y(i,:) = x(i) ./ z(i,:);
end
p = zeros(n,m);
sum_yj = sum(y);
for j = 1 : m
p(:,j) = y(:,j) / sum_yj(j);
end
e = zeros(m,1);
temp = 0.0;
for j = 1 : m
temp = 0.0;
for i = 1 : n
temp = temp + p(i,j) * log(p(i,j));
end
temp = - 1/(log(n)) * temp;
e(j) = temp;
end
w = zeros(m,1);
sum_ej = sum(1 - e);
for j = 1 : m
w(j) = (1 - e(j)) / sum_ej;
end
Z = [65.82 * 10^4 126.6 * 10^8 17.04 * 10^8 341 * 10^8 1.08 * 10^8 3 * 10^8]';
G = zeros(m,1);
temp = 0;
for j = 1 : m
for i = 1 : n
temp = temp + z(i,j) * (2 * sum(x(1:i)) - x(i)) / (Z(j) * 1.06 * 10^4);
end
G(j) = 1 - temp;
temp = 0.0;
end
temp = 0.0;
for j = 1 : m
temp = temp + w(j) * G(j);
end
f = temp;
function [c ceq] = mycon(x,z)
%
% 非线性约束条件
%
Z = [65.82 * 10^4 126.6 * 10^8 17.04 * 10^8 341 * 10^8 1.08 * 10^8 3 * 10^8]';
G0 = [0.4339 0.3781 0.3660 0.5003 0.3705 0.4445]';
m = 6;
n = 15;
G = zeros(m,1);
temp = 0;
for j = 1 : m
for i = 1 : n
temp = temp + z(i,j) * (2 * sum(x(1:i)) - x(i)) / (Z(j) * 0.018 * 10^4);
end
G(j) = 1 - temp;
temp = 0.0;
end
c = zeros(m,1);
for j = 1 : m
c(j) = G(j) - G0(j);
end
ceq = [];
%%
function [w G] = get_w_G(x,z)
%
%
[n m] = size(z);
y = zeros(n,m);
for i = 1 : n
y(i,:) = x(i) ./ z(i,:);
end
p = zeros(n,m);
sum_yj = sum(y);
for j = 1 : m
p(:,j) = y(:,j) / sum_yj(j);
end
e = zeros(m,1);
temp = 0.0;
for j = 1 : m
temp = 0.0;
for i = 1 : n
temp = temp + p(i,j) * log(p(i,j));
end
temp = - 1/(log(n)) * temp;
e(j) = temp;
end
w = zeros(m,1);
sum_ej = sum(1 - e);
for j = 1 : m
w(j) = (1 - e(j)) / sum_ej;
end
Z = [65.82 * 10^4 126.6 * 10^8 17.04 * 10^8 341 * 10^8 1.08 * 10^8 3 * 10^8]';
G = zeros(m,1);
temp = 0;
for j = 1 : m
for i = 1 : n
temp = temp + z(i,j) * (2 * sum(x(1:i)) - x(i)) / (Z(j) * 0.018 * 10^4);
end
G(j) = 1 - temp;
temp = 0.0;
end