It is always interesting to read GAMS models written by someone else. There are probably three things one can observe:
- A nice formulation or concept that is useful to learn about.
- A bad implementation: something that really should not be done that way.
- A piece of code that is correct and defensible, but I would write it differently. This includes things like style, layout, formatting, etc.
Data
Set i 'accounts' / lab, h1, h2, p1, p2 /;
Alias (i,j);
Table xb(i,j) 'original estimate' lab h1 h2 p1 p2 lab 15 3 130 80 h1 na h2 na p1 15 130 20 p2 25 40 55 ;
Parameter tb(i) 'original totals' / lab 220, (h1,h2) na, p1 190, p2 105 / tw(i) 'wights for totals' xw(i,j) 'weights for cells';
tw(i) = 1; xw(i,j) = 1$xb(i,j); tw(i)$(mapval(tb(i)) = mapval(na)) = 0; xw(i,j)$(mapval(xb(i,j)) = mapval(na)) = 0;
display tw, xw; |
tw(i)$(tb(i)=na) = 0;
xw(i,j)$(xb(i,j)=na) = 0;
There are two audiences when you write code: one is the machine (in this case the GAMS compiler) and the other, more importantly, are human beings. It is very important to write code that is as readable as possible. In the fragment above this is obviously violated by obscuring the meaning of these assignments.
Nerdy detail: the second assignment cannot be written as:ij(i,j) = xb(i,j) and xb(i,j)<>na;as the rhs resolves for na values to na and 0 which is na.
Model equations
Variable x(i,j) 'estimated cells' t(i) 'estimated totals' dev 'deviations';
Equation rbal(i) 'row balance' cbal(j) 'column balance' devsqr 'definition of square deviations';
rbal(i).. t(i) =e= sum(j$xb(i,j), x(i,j));
cbal(j).. t(j) =e= sum(i$xb(i,j), x(i,j));
devsqr.. dev =e= sum((i,j)$xw(i,j), xw(i,j)*sqr(xb(i,j) - x(i,j))/xb(i,j)) + sum(i$tw(i), tw(i)*sqr(tb(i) - t(i))/tb(i)); |
Model bal / all /;
x.l(i,j) = xb(i,j)$xw(i,j); t.l(i) = tb(i)$tw(i);
solve bal using nlp minimizing dev; |
Conclusions
References
- sambal.gms : Social Accounting Matrix Balancing Problem, https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_sambal.html
- Zenios, S.A., Drud, A. and Mulvey, J.M. (1989), Balancing large social accounting matrices with nonlinear network programming. Networks, 19: 569-585.
Appendix A: original sambal.gms
$title Social Accounting Matrix Balancing Problem (SAMBAL,SEQ=77)
$onText A Social Accounting Matrix captures all the circular flows in an economy and is called balanced if the row total equal the column totals. A sample problem illustrates the use of Nonlinear Programming to balance such matrices.
Zenios, S A, Drud, A S, and Mulvey, J, Balancing some large Social Accounting Matrices with Nonlinear Programming. Tech. rep., Department of Civil Engineering, Princeton University, 1986.
Keywords: nonlinear programming, social accounting matrix, statistics $offText
Set i 'accounts' / lab, h1, h2, p1, p2 /;
Alias (i,j);
Table xb(i,j) 'original estimate' lab h1 h2 p1 p2 lab 15 3 130 80 h1 na h2 na p1 15 130 20 p2 25 40 55 ;
Parameter tb(i) 'original totals' / lab 220, (h1,h2) na, p1 190, p2 105 / tw(i) 'wights for totals' xw(i,j) 'weights for cells';
tw(i) = 1; xw(i,j) = 1$xb(i,j); tw(i)$(mapval(tb(i)) = mapval(na)) = 0; xw(i,j)$(mapval(xb(i,j)) = mapval(na)) = 0;
display tw, xw;
Variable x(i,j) 'estimated cells' t(i) 'estimated totals' dev 'deviations';
Equation rbal(i) 'row balance' cbal(j) 'column balance' devsqr 'definition of square deviations';
rbal(i).. t(i) =e= sum(j$xb(i,j), x(i,j));
cbal(j).. t(j) =e= sum(i$xb(i,j), x(i,j));
devsqr.. dev =e= sum((i,j)$xw(i,j), xw(i,j)*sqr(xb(i,j) - x(i,j))/xb(i,j)) + sum(i$tw(i), tw(i)*sqr(tb(i) - t(i))/tb(i));
Model bal / all /;
x.l(i,j) = xb(i,j)$xw(i,j); t.l(i) = tb(i)$tw(i);
solve bal using nlp minimizing dev; |
Appendix B: proposed model
$onText
Matrix Balancing Problem Find a nearby matrix such that row and column totals are correct. Additional condition: preserve zeros, positive values. Slightly different relative quadratic objective
$offText
*--------------------------------------------------------------------- * data entry *---------------------------------------------------------------------
Sets iext 'accounts+totals' / lab, h1, h2, p1, p2, total / i(iext) 'accounts' / lab, h1, h2, p1, p2 / ;
Alias (iext,jext),(i,j);
Table xb(iext,jext) 'original estimate' lab h1 h2 p1 p2 total lab 15 3 130 80 220 h1 na na h2 na na p1 15 130 20 190 p2 25 40 55 105 total 220 na na 190 105 ;
* totals should be symmetric abort$sum(i$(xb(i,"total")<>xb("total",i)),1) "totals in XB not symmetric";
*--------------------------------------------------------------------- * derived data *---------------------------------------------------------------------
Parameter tb(i) 'totals'; tb(i) = xb(i,'total'); Sets ii(i) 'indicates if i is included in obj' ij(i,j) 'indicates if (i,j) is included in obj' nz(i,j) '(i,j) is nonzero: include in constraints' ;
ii(i) = tb(i)<>na; ij(i,j) = xb(i,j)<>0 and xb(i,j)<>na; nz(i,j) = xb(i,j); display ii, ij, nz;
* assumption: non-negative values abort$sum(ij$(xb(ij)<0),1) "negative values in XB"; abort$sum(ii$(tb(ii)<0),1) "negative values in TB";
*--------------------------------------------------------------------- * model *---------------------------------------------------------------------
Positive variable x(i,j) 'estimated cells' t(i) 'estimated totals'; Variable dev 'deviations'; Equation rbal(i) 'row balance' cbal(j) 'column balance' devsqr 'definition of squared deviations';
rbal(i).. t(i) =e= sum(nz(i,j), x(i,j));
cbal(j).. t(j) =e= sum(nz(i,j), x(i,j));
* this is a little bit different from the original model * I tend to use this form devsqr.. dev =e= sum(ij, sqr(x(ij)/xb(ij)-1)) + sum(ii, sqr(t(ii)/tb(ii)-1));
Model sambal /all/;
*--------------------------------------------------------------------- * solve *---------------------------------------------------------------------
* initial values x.l(ij) = xb(ij); t.l(ii) = tb(ii);
solve sambal using qcp minimizing dev;
*--------------------------------------------------------------------- * reporting *---------------------------------------------------------------------
parameter sol(iext,jext) 'final estimate'; sol(i,j) = x.l(i,j); sol(i,'total') = t.l(i); sol('total',i) = t.l(i); display xb,sol; |
No comments:
Post a Comment