% 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