Large scale Optimization using Sparse?

Hi everyone: I am working on large scale optimization using LINPROG, it has huge matrices of A,B,Aeq,Beq, it goes out of memory when I run it. The Matrices here have almost all zero elements in each row except for one or two columns in that row. I think using sparse matrix would help me resolve this problem but do I have to initialize the matrix as sparse before I start filling it up or after I have done so. Here matrices are filled as separate blocks so I am generalizing the row and column to be filled can this generalizing be done using sparse.? If anyone has a link or something illustrating an example for sparse matrix for large scale optimization please forward it to me at [e-mail]

Answer 1

[Gabo]

You need to create a sparse matrix directly. If you try to create a full matrix and then convert it to sparse, you'll run out of memory. It has happened to me with large matrices.

See the examples in sparse. You can specify the location of the nonzero entries ( i and j ), the nonzero values ( s ), the size of the matrix ( m and n ) and the maximum number of nonzero elements in it ( nzmax ).

[Samir]

Thanks..You are right I tried doing the sparse in the end it went out of memory before reaching the statement itself but when I do it initially I need to fill up the matrix using indexing. I have pasted the code below just paste on your matlab editor window it RED line on lines 62, 67, 101, 102(aeq matrix) when cursor kept on the highlighted aeq it shows sparse indexing likely to be slow.

clear;

clc;

%% Main system data

nb=3;

nl=3;

h=4;

%Storage and gen assumed at every bus

ng=nb;

nst=nb;

%_______________________________________________

%% Read line data

%Read entire table of line data

B=[1;1;2];

D=[2;3;3];

J=[0.01;0.02;0.03];

%__________________________________________________

%% Formation of Aeq

%aeq=zeros(nb*h+nl*h+1, ng*h+nl*h+nst*h+(nb-1)*h);

aeq=sparse(nb*h+nl*h+1, ng*h+nl*h+nst*h+(nb-1)*h);

%__________________________________________________

%Formation of bus dictionary

busdict=zeros(nb,1);

busdict(1)=1;

running=1;

for iline=1:nl;

ifrom = B(iline);

ito = D(iline);

xline=J(iline);

ifound=0;

for i=1:running;

if busdict(i)==ifrom;

ifound=i;

ifromn=i;

end;

end;

% if ifound is still zero, a new bus name was found and need to add

% this to the busdict

if ifound == 0;

running=running+1;

busdict(running)=ifrom;

ifromn=running;

end;

ifound=0;

for i=1:running;

if busdict(i)==ito;

ifound=i;

iton=i;

end;

end;

if ifound == 0

running = running+1;

busdict(running)=ito;

iton=running;

end;

%Form the part of aeq thaty deals with the line impedances -- processing

%one line at a time, working on line iline now

%No entry if starting from the reference bus

if ifromn ~= 1;

%add in -1/xline into aeq

%generate an h by h eye matrix multiplied by -1/x

blok=(-1/xline)*eye(h,h);

aeq(nb*h+ifromn-1+(iline-1)*h:nb*h+ifromn+h-2+(iline-1)*h, (nb+nst+nl)*h+1:(nb+nst+nl)*h+h)=blok;

end;

%No entry if ending at the reference bus

if iton~= 1;

blok=(1/xline)*eye(h,h);

aeq(nb*h+(iline-1)*h+1:nb*h+(iline-1)*h+h,(nb+nst+nl)*h+(iton-2)*h+1:(nb+nst+nl)*h+(iton-2)*h+h)=blok;

end;

end;

%________________________________________________________________________%

%Generator injection power - assumed at all busses. Zero generation

%accommodated at LB and UB -- continue on the formation of aeq

inb=eye(nb*h,nb*h);

aeq(1:nb*h,1:nb*h)=inb;

%__________________________________________________________________

%Storage power - assumed at all buses. Zero storage accomodated at LB n UB

aeq(1:nb*h,(nb)*h+1:(nb+nst)*h)=-inb;

%_____________________________________________________________________

%Stored Energy at the end of the day should be zero

aeq(nb*h+nl*h+1,(nb)*h+1:(nb+nst)*h)=24/h;

%______________________________________________________________________

%Line injection power

inh=eye(h,h);

for k=1:nl;

stb=B(k);

stbb=0;

for look=1:nb;

if busdict(look)==stb;

stbb=look;

end;

end;

endb=D(k);

endbb=0;

for look=1:nb;

if busdict(look)==endb;

endbb=look;

end;

end;

% %______________________________________________________________________

% %Line goes from stb to endb

aeq(h*(stbb-1)+1:h*(stbb-1)+h,nb*h+nst*h+1+h*(k-1):nb*h+nst*h+h*(k-1)+h)=-inh;

aeq(h*(endbb-1)+1:h*(endbb-1)+h,nb*h+nst*h+1+h*(k-1):nb*h+nst*h+h*(k-1)+h)=inh;

end;

%Power Flow vs. delta - this completes the formation of matrix aeq

inl=eye(nl*h,nl*h);

aeq(nb*h+1:(nb+nl)*h,nb*h+nst*h+1:(nb+nst+nl)*h)=inl;

Answer 2

[Gabo]

Without looking into all the details, the Code Analyzer seems to be recommending that you indicate the nonzero patter from the start, instead of declaring aeq as an empty sparse matrix, and then filling in the details.

Quick example. If your matrix looked like this:

aeq= [ 0.1   0   0   0   0;
        0   0.1  0   0   0;
        0    0   0  0.2  0;
        0    0   0   0  0.2]

you can do

i=[1 2 3 4];
j=[1 2 4 5];
s=[0.1 0.1 0.2 0.2];
m=4;
n=5;
aeq = sparse(i,j,s,m,n);

In other words, maybe you can use the for loops to create the vectors i,j,s, and create the sparse matrix after the loops. By the way, note that instead of filling in a square block with a full identity matrix, you can save memory by filling in the diagonal elements only.