$TITLE  QAP Example using penalty and circle constr, 5 by 5
$OFFUPPER
* This problem is the quadratic assignment problem, i.e. the
*  problem of locating facilities in sites given 
*  distance and flow matrices, as well as location costs
*
*  References: 
*       Burkard and Derigs book
*

  SETS
       I   sites / S1 * S5 /
       J   facilities   / F1 * F5 / ;

ALIAS (I,K);

ALIAS (J,L);


  TABLE A(I,I)  distance matrix between sites
    S1  S2  S3  S4  S5  
S1  7   1   1   2   3
S2  1   7   2   1   2
S3  1   2   7   1   2
S4  2   1   1   7   1
S5  3   2   2   1   7   ;

  TABLE B(J,J)  flows between facilities
   F1  F2  F3  F4  F5
F1 12  5   2   4   1
F2 5  12   3   0   2
F3 2   3  12   0   0
F4 4   0   0  12   5
F5 1   2   0   5  12 ;



  TABLE C(I,J)  location costs of plant j in location i 
     F1  F2  F3  F4  F5
S1  -84 -84 -84 -84 -84
S2  -84 -84 -84 -84 -84
S3  -84 -84 -84 -84 -84
S4  -84 -84 -84 -84 -84
S5  -84 -84 -84 -84 -84 ;


  VARIABLES
       X(I,J)  0 or 1 indicating facility j in location i
       PN      value of the constraint in the penalty function
       ZT      true objective value without penalties
       Z       total cost for objective function ;

  POSITIVE VARIABLE X ;
* X.L(I,J) = 1 ;

SCALAR P penalty parameter /1000/ ;

   X.L('S1','F2') = 1 ; 
   X.L('S2','F1') = 1 ; 
   X.L('S3','F3') = 1 ; 
   X.L('S4','F4') = 1 ; 
   X.L('S5','F5') = 1 ; 
 
* DISPLAY X.L ;

  EQUATIONS
       COST        define total objective function
       CONSTR1     single ball constraint  ;

  COST .. Z=E= SUM((I,J),C(I,J)*X(I,J)+SUM((K,L),A(I,K)*B(J,L)*X(I,J)*X(K,L))) 
               + P*(SUM( J,     SQR( SUM(I,X(I,J))  -1) ) 
                +   SUM( I,     SQR( SUM(J,X(I,J))  -1) )) ;
*               -(1/P)*SUM((I,J), LOG(X(I,J))  ) ;
            
  CONSTR1.. SUM(I,   SUM(J, SQR( X(I,J) )  ) ) =E= 5 ;

MODEL PNLTY2 /ALL/;

 SOLVE PNLTY2  MINIMIZING Z USING NLP;

  DISPLAY X.L, X.M ;

  ZT.L = SUM((I,J),C(I,J)*X.L(I,J)+SUM((K,L),A(I,K)*B(J,L)*X.L(I,J)*X.L(K,L))) ;

  PN.L = SUM(J,SQR(SUM(I,X.L(I,J))-1))+SUM(I,SQR(SUM(J,X.L(I,J))-1));

  DISPLAY ZT.L, PN.L ;