Introduction Ax b GAMS A - Amsterdam Optimization

1 SOLVING SYSTEMS OF LINEAR EQUATIONS WITH GAMSERWIN document describes some issues with respect to solving sys-tems of linear equationsAx=busing a system of linear equationsAx=b, whereAis a (square) matrix, is ataskGAMSand LP solvers are well equipped to do. As LP solvers do not requireAto be square, even more general systems can be solved. Even if the matrixAissingular, an LP solver will report a solution, provided there exists a trick that needs to be used to shoe horn such a model into a linear opti-mization model is to select a variable to optimize. One can just pick a variable thatappears in the model, especially when it is known that there exists just one uniquesolution. However, from a modeling point of view it may be better to introduce adummyvariable. This makes clear that this variable has no intrinsic meaning inthe model, but just is created for syntactical reasons. Below is a small example ofa system of linear equations that is not square.

2 We imposed bounds onpi, whichhelps to detect modeling errors: the model will be declared infeasible ifpi<0 orpi> $ontextGiven a Markov transition matrix, find the steady stateprobabilities p(i).Erwin Kalvelagen, 2001$offtextset i states /state1*state4/;alias (i,j);table A(i,j) transition matrix state1 state2 state3 state4state1 ;variablesdummy dummy objective variable p(i) steady state probabilities ;Date: February KALVELAGEN positive variables p; (i) = 1;equationssteadystate(i)normalizeedummy ;steadystate(i).. sum(j, A(i,j)*p(j)) =e= p(i); sum(i, p(i)) =e= 1; dummy =e= 0;model m1 /steadystate,normalize,edummy/;solve m1 minimizing dummy using lp;Note that if we have a square modelAx=bwith a nonsingular matrix we cansolve the system of equations very efficiently: namely in zero iterations. We knowthat all variables are basic and all equations are non-basic in the final an advanced basis, we can tell the solver to use such an optimal $ontextSolve3x + 2y + 4z = 54x - 3y - 7z = 172x - 3y + 6z = 2in zero iterations.

3 $offtextvariables x,y,z,dummy;equations e1,e2,e3,edummy; 3*x + 2*y + 4*z =e= 5; 4*x - 3*y - 7*z =e= 17; 2*x - 3*y + 6*z =e= 2; dummy =e= 0;** the variables are basic* ; ; ; ;** the equations are nonbasic* = 1; = 1; = 1; = 1;model m /all/;solve m using lp minimizing dummy;Indeed this model solves in zero iterations:2 O L V E S U M M A R YMODEL m OBJECTIVE dummyTYPE LP DIRECTION MINIMIZESOLVER BDMLP FROM LINE 40** SOLVER STATUS 1 NORMAL COMPLETION** MODEL STATUS 1 OPTIMAL** OBJECTIVE VALUE USAGE, LIMIT COUNT, LIMIT 0 10000If the matrixAis not square, or ifAis singular, the above approach can stillbe used, but the algorithm will need a few iterations to move some slack variablesinto the the inverseTo find the inverse of a matrix, we can use (or misuse) the looping facilities inGAMSto implement a Gaussian elimination procedure.

4 An example of aGAMS model that implements this strategy can be found in theGAMS model different approach is to find the inverse column by column. We can solve:(1)Axi=eiwhereeiis theithcolumn of the identity. Then the columns (x1,..,xn) form theinverse different approach is to write down(2)AX=Ias a set of linear equations that need to be solved forX=A that implements this, can be written as:Model $ontextFinds the inverse of a matrixErwin Kalvelagen, nov 2002 Reference: model from the model $offtextset i /i1*i3 /;alias (i,j,k);table a(i,j) original matrix i1 i2 i3i1 1 2 3i2 1 3 4i3 1 4 3;parameter ident(i,j) identity matrix ;ident(i,i) = 1;3 KALVELAGEN** method 1: solve a series of models Ax=b where* b is a column from the unit matrix.*variabledummy dummy objective variable col(j) column of inv(A) ;equationsedummy dummy objective function lineq(i) system of linear equations ;parameter b(i) rhs ; dummy =e= 0;lineq(i).

5 Sum(k, a(i,k)*col(k)) =e= b(i);model m1 /lineq,edummy/;parameter inverse(i,j) inv(A) ;loop(j,b(i) = ident(i,j);solve m1 minimizing dummy using lp;inverse(i,j) = (i););display inverse;** method 2: solve a bigger system AX=I*variable ainv(i,j) inverse of A ;equation lineqall(i,j);lineqall(i,j).. sum(k, a(i,k)*ainv(k,j)) =e= ident(i,j);model m2 /lineqall,edummy/;solve m2 minimizing dummy using lp;display ;The resulting solution for(3)A= 1 2 31 3 41 4 3 is given by:---- 42 VARIABLE inverse of Ai1 i2 i3i1 of complex linear , the complex all the solvers underneathGAMSare not equipped to deal with complex numbers. Complex numbers havethe formz=a+ibwithidefined by the famous expressioni2= called thereal part ofzandbforms the imaginary part. We can map a complex system ofequations inCnonto a problem inR2n. This allows us to solve the problem usingstandard system of equations with complex variables and coefficientsAx=bcan beexpanded into:(4)(re(A) im(A)im(A)re(A))(re(x)im(x))=(re(b)im(b) )As an example consider the problem:(1 +i)z+ (2 i)w= 2 + 7i7z+ (8 2i)w= 4 9i(5)The following model shows a simple implementation:Model $ontextSolution of a system of complex equationsErwin Kalvelagen, April 2002$offtextset i /i1*i2/;alias (i,j);table data(*,i,*)i1 i2 1 2 7 8 1 -1 -2 -9;variablesrx(i) real part ix(i) imaginary part dummy dummy objective variable ;equationsreal(i) real part imag(i) imaginary part dummyobj dummy objective ; dummy=e=0;real(i).

6 Sum(j, data( real ,i,j)*rx(j)-data( imag ,i,j)*ix(j)) =e= data( real ,i, rhs );imag(i)..sum(j, data( imag ,i,j)*rx(j)+data( real ,i,j)*ix(j)) =e= data( imag ,i, rhs );5 KALVELAGEN model m /all/;solve m using lp minimizing dummy;display , ;The solution of this system is:z=838185 699185iw= 698185+229185i(6)The solution of theGAMS model is identical to this solution:---- 47 VARIABLE real parti1 , i2 47 VARIABLE imaginary parti1 , i2 a complex same approach can be used to find the(complex) inverse of a complex matrix. We searchA 1such that:(7)(re(A) im(A)im(A)re(A))(re(A 1)im(A 1))=(I0)which can be written out as: kre(ai,k)re(a 1k,j) im(ai,k)im(a 1k,j) =ei,j kim(ai,k)re(a 1k,j) +re(ai,k)im(a 1k,j) = 0(8)The model below calculates the inverse matrix of(9)(1 +i2 i78 2i)Model $ontextUsing the example from ,the code below first calculates the complex inverse ainv,and then calculates x = inv(a)* Kalvelagen, May 2004$offtextset i /i1*i2/;alias (i,j,k);table data(*,i,*)i1 i2 1 2 7 8 1 -1 -2 -9;6 cmplx /re,im/;parameter a(i,j,cmplx);a(i,j, re ) = data( real ,i,j);a(i,j, im ) = data( imag ,i,j);parameter ident(i,j) identity matrix ;ident(i,i) = 1;parameter b(i,cmplx);b(i, re ) = data( real ,i, rhs );b(i, im ) = data( imag ,i, rhs );display a,ident;variablesainv(i,j,cmplx) inverse matrix dummy dummy objective variable ;equationsreal(i,j)imag(i,j)dummyobj dummy objective ; dummy=e=0;real(i,j).

7 Sum(k, a(i,k, re )*ainv(k,j, re )-a(i,k, im )*ainv(k,j, im )) =e= ident(i,j);imag(i,j)..sum(k, a(i,k, im )*ainv(k,j, re )+a(i,k, re )*ainv(k,j, im )) =e= 0;model m /all/;solve m using lp minimizing dummy;display ;parameter x(i,cmplx);x(i, re ) = sum(j, (i,j, re )*b(j, re ) - (i,j, im )*b(j, im ));x(i, im ) = sum(j, (i,j, re )*b(j, im ) + (i,j, im )*b(j, re ));display b,x; s a of Partial Differential Equations (PDE s) us-ing finite difference approximation scheme s leads to large systems of linear equa-tion. Consider the numerical solution of the Laplace s Equation:(10) 2u x2+ 2u y2= 0on a square, where the boundary values ofuare given. Assuming a grid of theform:(11)ui,j+1ui 1,jui,jui+1,jui,j 1then the first derivatives can be approximated by: u x ui+1,j ui,jh u y ui,j+1 ui,jh(12)8 ERWIN KALVELAGENand the second derivatives by: 2u x2 1h(ui+1,j ui,jh ui,j ui 1,jh)=ui+1,j 2ui,j+ui 1,jh2 2u y2 1h(ui,j+1 ui,jh ui,j ui,j 1h)=ui,j+1 2ui,j+ui,j 1h2(13)Combining this we can approximate(14) 2u x2+ 2u y2= 0by(15)1h2(ui,j+1+ui,j 1+ui 1,j+ui+1,j 4ui,j) = 0or(16)4ui,j=ui,j+1+ui,j 1+ui 1,j+ui+1,jAs an example consider the problem of solving the steady-state heat conductionproblem over the unit square.

8 This can be formulated as solving Laplace s Equationsubject to the Dirichlet boundary conditions:u(x,0) = 300u(x,1) =u(0,y) =u(1,y) = 0(17)This model can be implemented inGAMSas follows:Model $ontextFinite difference method of Laplace s EquationUxx + Uyy = 0subject to the Dirichlet boundary conditionErwin Kalvelagen, april 2003$offtextscalar b boundary value for y=0 /300/;set i interior nodes /i1*i99/;alias (i,j);set j0(j);j0(j)$(ord(j)=1) = yes;variables u(i,j),dummy;equation e(i,j),obj; dummy =e= 0;e(i,j).. 4*u(i,j) =e= u(i,j+1) + u(i,j-1) + u(i-1,j) + u(i+1,j) + b$(j0(j));model m /all/;7 m using lp minimizing dummy;display ; is quite a difficult LP, and many solvers will encounternumerical difficulties. The solution is again to provide a good initial ; (i,j)=1; ; (i,j)=0;This will significantly speed up the solution times. for a few solvers we havethe following results:SolverDefaultAdvanced solvers really appreciate the advanced basis resulting in a significant im-provement in 3D plot of the results the can export the solution to a GNUPLOT data fileusing:parameter v(i);v(i) = ord(i)/(1+card(i));file f ;put f;loop((i,j),put v(i), v(j), (i,j)/;);The result is shown in figure large input-output analysis deals with the inter-relationshipsbetween sectors of an economy.

9 Sectors use intermediate products from other prod-ucts and deliver goods to other sectors. Monetary flows go the other way I-O table records the transactions between sectors. Often they are organizedon yearly data, and several aggregation levels. A small I-O table deals with a hand-ful of sectors while a highly disaggregated I-O table may record the transactions ofhundreds of basic structure of an input-output transaction table is depicted in table ,1 x1, ,1 xn, VnDomesticProductionX1 XnTable table structureThe entriesxi,jcan be interpreted as the output of sectorithat is used by sectorjas input. The following relations hold:(18) jxi,j+Yi=Xiand(19) ixi,j+Vj=XjThe matrix of input-coefficientsAis defined by(20)ai,j=xi,j/XjThis leads to the equation(21)AX+Y=Xor(22)(I A)X=Y11which gives:(23)X= (I A) 1 YThe matrix (I A) 1is often called theLeontief Inverse Matrix. Calculatingthis inverse can be accomplished inGAMSby defining a system of linear equations(24)(I A) 1I= (I A)A different way is to use the expansion:(25)(I A) 1=I+A+A2+A3 We can easily derive this series expansion by writing:(I A)(I+A+A2+A3 ) =(I+A+A2+A3 ) (A+A2+A3 ) =I(26)if limk Ak= series converges as we have elements with a range ( 1,1).

10 In fact for thelarge example below, the maximum element ofAkis as follows:441 PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = PARAMETER maxelem = the series converges matrix (I A) 1is convenient to evaluate scenario s(27) X= (I A) 1 real-world example: Input-Output table for Japan, illus-trate the above we downloaded a real-world 93 sector I-O table from Japan for1995 [3]. The data came in form of an Excel spreadsheet, so we usedXLS2 GMSto convert this to aGAMS include file. The include file is too large to completelyreproduce here, but a fragment is listed below:* ---------------------------------------- -------------* XLS2 GMS Version , December 2001* Erwin Kalvelagen, GAMS Development Corp.