Skip to content

Commit b931eda

Browse files
siyabcsiyabc
authored
Add files via upload
1 parent 451932d commit b931eda

File tree

4 files changed

+330
-0
lines changed

4 files changed

+330
-0
lines changed

‎con2vert.m

Lines changed: 183 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
function [V,nr] = con2vert(A,b)
2+
% CON2VERT - convert a convex set of constraint inequalities into the set
3+
% of vertices at the intersections of those inequalities;i.e.,
4+
% solve the "vertex enumeration" problem. Additionally,
5+
% identify redundant entries in the list of inequalities.
6+
%
7+
% V = con2vert(A,b)
8+
% [V,nr] = con2vert(A,b)
9+
%
10+
% Converts the polytope (convex polygon, polyhedron, etc.) defined by the
11+
% system of inequalities A*x <= b into a list of vertices V. Each ROW
12+
% of V is a vertex. For n variables:
13+
% A = m x n matrix, where m >= n (m constraints, n variables)
14+
% b = m x 1 vector (m constraints)
15+
% V = p x n matrix (p vertices, n variables)
16+
% nr = list of the rows in A which are NOT redundant constraints
17+
%
18+
% NOTES: (1) This program employs a primal-dual polytope method.
19+
% (2) In dimensions higher than 2, redundant vertices can
20+
% appear using this method. This program detects redundancies
21+
% at up to 6 digits of precision, then returns the
22+
% unique vertices.
23+
% (3) Non-bounding constraints give erroneous results; therefore,
24+
% the program detects non-bounding constraints and returns
25+
% an error. You may wish to implement large "box" constraints
26+
% on your variables if you need to induce bounding. For example,
27+
% if x is a person's height in feet, the box constraint
28+
% -1 <= x <= 1000 would be a reasonable choice to induce
29+
% boundedness, since no possible solution for x would be
30+
% prohibited by the bounding box.
31+
% (4) This program requires that the feasible region have some
32+
% finite extent in all dimensions. For example, the feasible
33+
% region cannot be a line segment in 2-D space, or a plane
34+
% in 3-D space.
35+
% (5) At least two dimensions are required.
36+
% (6) See companion function VERT2CON.
37+
% (7) ver 1.0: initial version, June 2005
38+
% (8) ver 1.1: enhanced redundancy checks, July 2005
39+
% (9) Written by Michael Kleder
40+
%
41+
% EXAMPLES:
42+
%
43+
% % FIXED CONSTRAINTS:
44+
% A=[ 0 2; 2 0; 0.5 -0.5; -0.5 -0.5; -1 0];
45+
% b=[4 4 0.5 -0.5 0]';
46+
% figure('renderer','zbuffer')
47+
% hold on
48+
% [x,y]=ndgrid(-3:.01:5);
49+
% p=[x(:) y(:)]';
50+
% p=(A*p <= repmat(b,[1 length(p)]));
51+
% p = double(all(p));
52+
% p=reshape(p,size(x));
53+
% h=pcolor(x,y,p);
54+
% set(h,'edgecolor','none')
55+
% set(h,'zdata',get(h,'zdata')-1) % keep in back
56+
% axis equal
57+
% V=con2vert(A,b);
58+
% plot(V(:,1),V(:,2),'y.')
59+
%
60+
% % RANDOM CONSTRAINTS:
61+
% A=rand(30,2)*2-1;
62+
% b=ones(30,1);
63+
% figure('renderer','zbuffer')
64+
% hold on
65+
% [x,y]=ndgrid(-3:.01:3);
66+
% p=[x(:) y(:)]';
67+
% p=(A*p <= repmat(b,[1 length(p)]));
68+
% p = double(all(p));
69+
% p=reshape(p,size(x));
70+
% h=pcolor(x,y,p);
71+
% set(h,'edgecolor','none')
72+
% set(h,'zdata',get(h,'zdata')-1) % keep in back
73+
% axis equal
74+
% set(gca,'color','none')
75+
% V=con2vert(A,b);
76+
% plot(V(:,1),V(:,2),'y.')
77+
%
78+
% % HIGHER DIMENSIONS:
79+
% A=rand(15,5)*1000-500;
80+
% b=rand(15,1)*100;
81+
% V=con2vert(A,b)
82+
%
83+
% % NON-BOUNDING CONSTRAINTS (ERROR):
84+
% A=[0 1;1 0;1 1];
85+
% b=[1 1 1]';
86+
% figure('renderer','zbuffer')
87+
% hold on
88+
% [x,y]=ndgrid(-3:.01:3);
89+
% p=[x(:) y(:)]';
90+
% p=(A*p <= repmat(b,[1 length(p)]));
91+
% p = double(all(p));
92+
% p=reshape(p,size(x));
93+
% h=pcolor(x,y,p);
94+
% set(h,'edgecolor','none')
95+
% set(h,'zdata',get(h,'zdata')-1) % keep in back
96+
% axis equal
97+
% set(gca,'color','none')
98+
% V=con2vert(A,b); % should return error
99+
%
100+
% % NON-BOUNDING CONSTRAINTS WITH BOUNDING BOX ADDED:
101+
% A=[0 1;1 0;1 1];
102+
% b=[1 1 1]';
103+
% A=[A;0 -1;0 1;-1 0;1 0];
104+
% b=[b;4;1000;4;1000]; % bound variables within (-1,1000)
105+
% figure('renderer','zbuffer')
106+
% hold on
107+
% [x,y]=ndgrid(-3:.01:3);
108+
% p=[x(:) y(:)]';
109+
% p=(A*p <= repmat(b,[1 length(p)]));
110+
% p = double(all(p));
111+
% p=reshape(p,size(x));
112+
% h=pcolor(x,y,p);
113+
% set(h,'edgecolor','none')
114+
% set(h,'zdata',get(h,'zdata')-1) % keep in back
115+
% axis equal
116+
% set(gca,'color','none')
117+
% V=con2vert(A,b);
118+
% plot(V(:,1),V(:,2),'y.','markersize',20)
119+
%
120+
% % JUST FOR FUN:
121+
% A=rand(80,3)*2-1;
122+
% n=sqrt(sum(A.^2,2));
123+
% A=A./repmat(n,[1 size(A,2)]);
124+
% b=ones(80,1);
125+
% V=con2vert(A,b);
126+
% k=convhulln(V);
127+
% figure
128+
% hold on
129+
% for i=1:length(k)
130+
% patch(V(k(i,:),1),V(k(i,:),2),V(k(i,:),3),'w','edgecolor','none')
131+
% end
132+
% axis equal
133+
% axis vis3d
134+
% axis off
135+
% h=camlight(0,90);
136+
% h(2)=camlight(0,-17);
137+
% h(3)=camlight(107,-17);
138+
% h(4)=camlight(214,-17);
139+
% set(h(1),'color',[1 0 0]);
140+
% set(h(2),'color',[0 1 0]);
141+
% set(h(3),'color',[0 0 1]);
142+
% set(h(4),'color',[1 1 0]);
143+
% material metal
144+
% for x=0:5:720
145+
% view(x,0)
146+
% drawnow
147+
% end
148+
149+
c = A\b;
150+
if ~all(A*c < b);
151+
[c,f,ef] = fminsearch(@obj,c,'params',{A,b});
152+
if ef ~= 1
153+
error('Unable to locate a point within the interior of a feasible region.')
154+
end
155+
end
156+
b = b - A*c;
157+
D = A ./ repmat(b,[1 size(A,2)]);
158+
[k,v2] = convhulln([D;zeros(1,size(D,2))]);
159+
[k,v1] = convhulln(D);
160+
if v2 > v1
161+
error('Non-bounding constraints detected. (Consider box constraints on variables.)')
162+
end
163+
nr = unique(k(:));
164+
G = zeros(size(k,1),size(D,2));
165+
for ix = 1:size(k,1)
166+
F = D(k(ix,:),:);
167+
G(ix,:)=F\ones(size(F,1),1);
168+
end
169+
V = G + repmat(c',[size(G,1),1]);
170+
[null,I]=unique(num2str(V,6),'rows');
171+
V=V(I,:);
172+
return
173+
function d = obj(c,params)
174+
A=params{1};
175+
b=params{2};
176+
d = A*c-b;
177+
k=(d>=-1e-15);
178+
d(k)=d(k)+1;
179+
d = max([0;d]);
180+
return
181+
182+
183+

‎maxmin.m

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
2+
function[p,power_evolution]=maxmin(G,n,beta,pmax)
3+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4+
% Problem parameter initialization
5+
% L: number of users.
6+
% G: channel gain.
7+
% F: nonnegative matrix. F_{lj} = G_{lj} if l ~= j, and F_{lj} = 0 if l = j
8+
% v: nonnegative vector. v_l = 1/G_{ll}
9+
% beta: a vector that is a weight assigned to links to reflect priority.
10+
% pmax: upper bound of the total power constraints.
11+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
12+
F = diag(diag(1./G))*(G-diag(diag(G)));
13+
v = n./diag(G);
14+
[L,~] = size(G);
15+
p = rand(L,1);
16+
power_evolution = [p'];
17+
tol = 10e-7;
18+
err = 1;
19+
20+
while err > tol
21+
p_temp = p;
22+
sinr = p./(F*p+v);
23+
p = beta./sinr.*p;
24+
power_evolution = [power_evolution;(p/sum(p)*pmax)'];
25+
err = norm(p-p_temp,2);
26+
end
27+
p = p/sum(p)*pmax;
28+
end
29+

‎outage.m

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
function [p,power_evolution]=outage(G,n,beta,pmax)
2+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3+
% Problem parameter initialization
4+
% L: number of users.
5+
% G: channel gain.
6+
% F: nonnegative matrix. F_{lj} = G_{lj} if l ~= j, and F_{lj} = 0 if l = j
7+
% n: noise vector
8+
% v: nonnegative vector. v_l = n_l/G_{ll}
9+
% beta: a vector that is a weight assigned to links to reflect priority.
10+
% pmax: upper bound of the total power constraints.
11+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
12+
13+
F = diag(diag(1./G))*(G-diag(diag(G)));
14+
v = n./diag(G);
15+
[L,~] = size(G);
16+
p = rand(L,1);
17+
power_evolution = [p'];
18+
tol = 10e-7;
19+
err = 1;
20+
21+
while err > tol
22+
p_temp = p;
23+
p = (v.*beta./p + log(prod(1+diag(beta./p)*F*diag(p),2))).*p;
24+
power_evolution = [power_evolution;(p/sum(p)*pmax)'];
25+
err = norm(p-p_temp,2);
26+
end
27+
p = p/sum(p)*pmax;
28+
end
29+

‎outer_apprx.m

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
function [k,power,power_evolution]=outer_apprx(G,n,w,a,pmax)
2+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3+
% max_iteration: maximal number of iterations.
4+
% epsilon: stopping criterion if ||p(k+1)-p(k)|| <= epsilon.
5+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6+
F = diag(diag(1./G))*(G-diag(diag(G)));
7+
v = n./diag(G);
8+
max_iteration = 100;
9+
epsilon = 0.01;
10+
[L,~] = size(G);
11+
12+
% Initialize the polyhedral convex set.
13+
% A*gamma_tilde <= b is the polyhedral convex set.
14+
A = []; b = [];
15+
for l = 1:1:L
16+
B(:,:,l) = F+(1/pmax(l)).*(v*a(:,l)');
17+
18+
[vec lamb] = eig(B(:,:,l));
19+
[val index]=max(max(lamb));
20+
x = vec(:,index);
21+
[vec lamb] = eig(B(:,:,l).');
22+
[val index]=max(max(lamb));
23+
y = vec(:,index);
24+
xy = (x.*y./sum(x.*y));
25+
A = [A;xy'];
26+
b = [b;-log(max(eig(B(:,:,l))))];
27+
end
28+
A = [A; -eye(L)]; b=[b; 100*ones(L,1)];
29+
30+
gamma_tilde = rand(L,1);
31+
tolerance = 1;
32+
k = 0;
33+
power_evolution = [];
34+
35+
while tolerance >= epsilon || k <= max_iteration
36+
37+
gamma_tilde_k = gamma_tilde;
38+
39+
% Step 1 in Algorithm 1: compute the vertices.
40+
[V,nr] = con2vert(A,b);
41+
42+
% Step 2 in Algorithm 1: select max{sumrate} from the vertices.
43+
sumrate = -Inf;
44+
s = size(V,1);
45+
for t = 1:1:s
46+
temp = sum(w.*log(ones(L,1)+exp(V(t,:)')));
47+
if temp > sumrate
48+
sumrate = temp;
49+
gamma_tilde = V(t,:)';
50+
end
51+
end
52+
53+
% Step 3 in Algorithm 1: compute power.
54+
power = inv( eye(L) - diag(exp(gamma_tilde)) * F ) * diag(exp(gamma_tilde)) * v;
55+
power_evolution = [power_evolution ; power'];
56+
57+
% Step 4 in Algorithm 1: let J = max{ diag(exp(gamma_tilde))*B(:,:,l) }.
58+
J=1;
59+
rho = 0;
60+
for l = 1:1:L
61+
Btemp = diag(exp(gamma_tilde))*B(:,:,l);
62+
if max(eig(Btemp)) > rho
63+
rho = max(eig(Btemp));
64+
J = l;
65+
end
66+
end
67+
68+
% Step 5 in Algorithm 1: compute the Perron right and left eigenvectors
69+
% of diag(exp(gamma_tilde))*B(:,:,J).
70+
Btemp = diag(exp(gamma_tilde))*B(:,:,J);
71+
[vec lamb] = eig(Btemp);
72+
[val index]=max(max(lamb));
73+
x = vec(:,index);
74+
[vec lamb] = eig(Btemp.');
75+
[val index]=max(max(lamb));
76+
y = vec(:,index);
77+
xy = x.*y./sum(x.*y);
78+
% Step 6 in Algorithm 1: add a hyperplane to the polyhedral convex set.
79+
A = [A;xy'];
80+
b = [b;xy'*gamma_tilde-log(max(eig(Btemp)))];
81+
82+
% Step 7 in Algorithm 1: update the iteration number.
83+
tolerance = norm(gamma_tilde_k-gamma_tilde);
84+
k = k+1;
85+
end
86+
87+
% plot the power evolution.
88+
89+
end

0 commit comments

Comments
 (0)