% moqpgen: generate s p a r s e MOQP and write to ampl file.
%
% matlab by Sven Leyffer, Argonne National Laboratory, 2005.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Randomly generate an instance of an MOQP (multi-objective QP)
% and write the resulting model out to an AMPL file.
%
% The problem has p objective functions, n variables, and m constraints:
%
% minimize { x'*G_k*x/2 + g_k'*x , for k=1,...,p }
% subject to A*x >= b
% x >= 0 (not generated)
%
% where
% g_k, G_k objective (quadratic, pos. def.; size is n x n)
% A is m by n random
% b is m by 1 random >= 0 such that A*e >= b
%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 3; % ... number of objectives
m = 20; % ... number of constraints
n = 40; % ... number of variables
range_A = [-4, 4]; % ... map [0,1] into [-4,4] here
range_b = [ 0, 8]; % ... range of values for b
range_g = [-6, 6]; % ... range of values for g
dens_G = 0.02; % ... density of G_k
dens_A = 0.2; % ... density of A
outprec = 0; % ... output precision (0=2 digits, 1=double)
plotpics = 0; % ... 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
formatG = ['\nlet G[ %-6i , %-6i , %-6i ] := \t',outlenr,';'];
formatg = ['\nlet g[ %-6i , %-6i ] := \t',outlenr,';'];
formatA = ['\nlet A[ %-6i , %-6i ] := \t',outlenr,';'];
formatb = ['\nlet b[ %-6i ] := \t',outlenr,';'];
% ... some useful constants
e = ones(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATE S P A R S E DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ... generate objective Hessian [NB: BLOCK DIAGONAL]
G = zeros(n,n,p);
for i=1:p
% ... generate symmetric random matrix
G(:,:,i) = sprandsym(n,dens_G);
% ... make sure Hessians are positive definite
mineigG = eigs(G(:,:,i),1,'SA');
G(:,:,i) = G(:,:,i) + max(-mineigG,0)*speye(n,n);
% ... 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,p);
jG = zeros(nnzGmax,p);
vG = zeros(nnzGmax,p);
for i=1:p
[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(p,n);
for i=1:p
g(i,:) = range_g(1) + rand(n,1) * (range_g(2) - range_g(1));
end % for
% ... generate constraint matrices in given range & ensure that A*e >= b
A = zeros(m,n);
b = zeros(m,1);
AA = sprand(m,n,dens_A);
[iA,jA,vA] = find(AA);
vA = range_A(1) + vA * (range_A(2) - range_A(1));
A = full( sparse(iA,jA,vA,m,n) );
nnzA = nnz(A);
b = range_b(1) + rand(m,1) * (range_b(2) - range_b(1));
b = min( b , A*e );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISPLAY SPARSITY PATTERNS GENERATED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if plotpics
figure(1); clf; spy(A); title('sparsity pattern of A');
for i=1:p
figure(i+1); clf; spy(G(:,:,i)); title(['sparsity pattern of G ',num2str(i)]);
end % for
end % if
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT PROBLEM TO AMPL FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ... open file
fid = fopen('moqp.dat','w+');
fprintf (fid,'# random MOQP s p a r s e data from matlab\n');
% ... write dimensions of problem
fprintf(fid,'\nparam n := %-4i ; ',n); fprintf(fid,'# ... number of variables');
fprintf(fid,'\nparam m := %-4i ; ',m); fprintf(fid,'# ... number of constraints');
fprintf(fid,'\nparam p := %-4i ; ',p); fprintf(fid,'# ... number of objectives');
fprintf(fid,'\n');
fprintf(fid,'# Hessian matrices G\n');
for i=1:p
for j=1:nnzG(i)
fprintf(fid,formatG,i,iG(j,i),jG(j,i),vG(j,i));
end % for
fprintf(fid,' \n');
end % for
fprintf(fid,'# Linear objectives g\n');
for i=1:p
for j=1:n
fprintf(fid,formatg,i,j,g(i,j));
end % for
fprintf(fid,' \n');
end % for
fprintf(fid,'# Linear constraints matrix A\n');
for j=1:nnzA
fprintf(fid,formatA,iA(j),jA(j),vA(j));
end % for
fprintf(fid,'\n');
fprintf(fid,'# Linear constraints: rhs b\n');
for j=1:m
fprintf(fid,formatb, j, b(j));
end % for
fprintf(fid,'\n');
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% E N D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%