% epecgen: generate s p a r s e EPEC and write to ampl file.
%
% matlab by Sven Leyffer, Argonne National Laboratory, 2004.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Randomly generate an instance of an EPEC (equilibrium problem
% with equilibrium constraints) and write the resulting model out
% to an AMPL file.
%
% This is an instance of a multi-leader-follower game, with L leaders.
%
% The problem for leader l=1,...,L is (all have x_l \in R^n):
%
%
% minimize x'*G_l*x/2 + g_i'*x
% subject to A_l*x_l <= b_l
% 0 <= y _|_ N*x + M*y + q >= 0
%
% where leader l's problem is defined by
% g_l, G_l objective (quadratic, pos. def.; size is n x L)
% A_l is p by n random
% b_l is p by 1 random >= 0 such that A_l*e <= b_l
%
% and the followers' problem is defined by
% N is m by (n x L) random
% M is m by m random \in [0,1] diagonally dominant
% q is m by 1 random >= 0
%
% using a uniform (0,1) random distribution (poss. scaled).
% To use other distributions replace "rand(" by "randn(" etc.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPECIFY PROBLEM DIMENSIONS, DATA RANGE AND DENSITY HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l = 4; % ... number of leaders
p = 8; % ... number of control constraints
m = 32; % ... number of states y
n = 16; % ... number of controls x
range_A = [-4, 4]; % ... map [0,1] into [-4,4] here
range_N = [ 0, 8]; % ... map [0,1] into [ 0,8] here
range_b = [ 0, 8]; % ... range of values for b
range_g = [-6, 6]; % ... range of values for g
range_q = [-4, 4]; % ... range of values for q
dens_G = 0.02; % ... density of G_i
dens_A = 0.2; % ... density of A
dens_N = 0.05; % ... density of N
dens_M = 0.1; % ... density of M
outprec = 0; % ... output precision (0=2 digits, 1=double)
plotpics = 1; % ... provide plots of sparsity pattern or not
%rand('seed',0); % ... uncomment for different random seed
% ... precision of output formats for printed outputs
if (outprec)
outlenr = '%20.16g '; % ... double precision (16 digits)
else
outlenr = '%10.2g '; % ... low precision (2 digits)
end % if
form1i1g = ['\n\t%-6i \t',outlenr];
form2i1g = ['\n\t%-6i %-6i \t',outlenr];
form3i1g = ['\n\t%-6i %-6i %-6i \t',outlenr];
% ... some useful constants
nl = n*l;
e = ones(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATE S P A R S E DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ... generate objective Hessian [NB: BLOCK DIAGONAL]
G = zeros(nl,nl,l);
for i=1:l
% ... generate symmetric random matrix
G(:,:,i) = sprandsym(nl,dens_G);
% ... make sure Hessians are positive definite
mineigG = eigs(G(:,:,i),1,'SA');
G(:,:,i) = G(:,:,i) + max(-mineigG,0)*speye(nl,nl);
% ... find number of nonzeros for plots
nnzG(i) = nnz( G(:,:,i) );
end % for
nnzGmax = max( nnzG ); % ... find sparsity pattern of G's
iG = zeros(nnzGmax,l);
jG = zeros(nnzGmax,l);
vG = zeros(nnzGmax,l);
for i=1:l
[iG(1:nnzG(i),i),jG(1:nnzG(i),i),vG(1:nnzG(i),i)] = find(G(:,:,i));
end % for
% ... generate objective gradient in given range
g = zeros(nl,l);
for i=1:l
g(:,i) = range_g(1) + rand(nl,1) * (range_g(2) - range_g(1));
end % for
% ... generate leaders' constraint matrices in given range
A = zeros(p,n,l);
b = zeros(p,l);
for i=1:l
AA = sprand(p,n,dens_A);
[iA,jA,vA] = find(AA);
vA = range_A(1) + vA * (range_A(2) - range_A(1));
A(:,:,i) = full( sparse(iA,jA,vA,p,n) );
nnzA(i) = nnz(A(:,:,i));
% ... ensure that A*e <= b
b = range_b(1) + rand(p,1) * (range_b(2) - range_b(1));
Ae = A(:,:,i)*e;
b(:,i) = max( b , Ae );
end % for
nnzAmax = max( nnzA ); % ... find sparsity pattern of A's
iA = zeros(nnzAmax,l);
jA = zeros(nnzAmax,l);
vA = zeros(nnzAmax,l);
for i=1:l
[iA(1:nnzA(i),i),jA(1:nnzA(i),i),vA(1:nnzA(i),i)] = find(A(:,:,i));
end % for
% ... generate follower's LCP: matrix N
N = sprand(m,nl,dens_N);
[iN,jN] = find(N);
vN = range_N(1) + N * (range_N(2) - range_N(1));
nnzN = nnz(N);
% ... generate follower's LCP: matrix M & ensure it is diagonally dominant
M = sprand(m,m,dens_M);
dM = diag(M); % subtract diagonal from M
MM = sparse([1:m],[1:m],dM);
Md = M - MM;
col_sum = sum(Md,1); % work out column sum
row_sum = sum(Md,2); % work out row sum
max_sum = max( col_sum , row_sum' ) + 1E-2; % max of two above + 1E-2
for i=1:m
M(i,i) = max_sum(i);
end % for
nnzM = nnz(M);
[iM,jM,vM] = find(M);
% ... generate follower's LCP: right hand side q
q = range_q(1) + rand(m,1) * (range_q(2) - range_q(1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISPLAY SPARSITY PATTERNS GENERATED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if plotpics
figure(1); clf; spy(N); title('sparsity pattern of N');
figure(2); clf; spy(M); title('sparsity pattern of M');
for i=1:l
figure(2+2*i-1); clf; spy(A(:,:,i)); title(['sparsity pattern of A ',num2str(i)]);
figure(2+2*i); clf; spy(G(:,:,i)); title(['sparsity pattern of G ',num2str(i)]);
end % for
end % if
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT PROBLEM TO AMPL FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ... open file
fid = fopen('epec.dat','w+');
fprintf (fid,'# random EPEC s p a r s e data from matlab\n');
% ... write dimensions of problem
fprintf(fid,'\nparam l := %-4i ; ',l); fprintf(fid,'# ... number of leaders');
fprintf(fid,'\nparam n := %-4i ; ',n); fprintf(fid,'# ... number of leaders variables');
fprintf(fid,'\nparam m := %-4i ; ',m); fprintf(fid,'# ... number of followers');
fprintf(fid,'\nparam p := %-4i ; ',p); fprintf(fid,'# ... number of leaders constraints');
fprintf(fid,'\n');
% ... Hessian matrices of leaders
fprintf(fid,'\nparam:\t \tG :=');
for i=1:l
for j=1:nnzG(i)
fprintf(fid,form3i1g,iG(j,i),jG(j,i),i,vG(j,i));
end % for
if (i