$*************************************************************
$ Multi-commodity capacity facility location-allocation
$
$ Optimal Solution:  3.31837e+08
$*************************************************************

OPTION {{
  FTOL = "1e-6";
}}

DECLARATION {{
  INDEX {i,j,k,p};
  SET I = |1:3|;
  SET J = |1:3|;
  SET K = |1:3|;
  SET P = |1:2|;

  PARAM c(I,J,K,P) = {115.274 ,  98.5559, 142.777 ,  33.9886, 163.087 ,  10.4376,
                      234.406 , 142.066 ,  50.6436, 123.61  , 242.356 , 135.071 ,
                       10.7347,  56.0272,  14.912 , 169.218 , 209.028 , 259.29  ,
                      165.41  ,  40.7497, 124.907 ,  18.495 ,  95.2789, 251.899 ,
                      114.185 ,  37.8148,  10.5547,  52.5162,  37.4727, 254.843 ,
                      266.645 , 136.583 ,  15.092 , 194.101 ,  78.768 , 120.36  ,
                      257.318 , 172.747 , 142.813 , 251.331 ,  15.9113,  48.8251,
                      289.116 , 129.705 , 275.621 ,  20.2235, 253.789 ,  56.7474,
                      201.646 , 164.573 , 295.157 , 151.474 , 221.794 , 278.304 };

  PARAM s(I,P) = {60, 60, 60, 60, 60, 60};
  PARAM d(K,P) = {60, 60, 60, 60, 60, 60};
  PARAM f(J) = {2.4814e+06, 2.15646e+06, 2.09773e+06};
  PARAM v(J) = {276.28, 792.912, 991.679};
  PARAM vu(J) = {2749.5, 2872.94, 2508.06};
  PARAM vl(J) = {50, 50, 50};

  XVAR {x(I,J,K,P)};
  YVAR {y(J,K),z(J)};
  BINARY {y(J,K),z(J)};
  POSITIVE {x(I,J,K,P)};
  UBDS x(I,J,K,P) = <i E I| <j E J| <k E K| <p E P| 1000> > > >;
  YSTP {1, 0, 0,
        0, 1, 0,
        0, 0, 1,
        1, 1, 1};
}}

MODEL {{
  MIN: <<i E I| <<j E J| <<k E K| <<p E P| c(i,j,k,p)*x(i,j,k,p) >> >> >> >>
       + <<j E J| f(j)*z(j) >>
       + <<j E J| v(j)*<<i E I| <<k E K| <<p E P| x(i,j,k,p) >> >> >>^2.5>>;
  lsum1(i E I, p E P): <<j E J| <<k E K| x(i,j,k,p) >> >> =l= s(i,p);
  lsum2(j E J, k E K, p E P): <<i E I| x(i,j,k,p) >> =E= d(k,p)*y(j,k);
  lsum3(j E J): <<k E K| <<p E P| d(k,p)*y(j,k) >> >> =l= vu(j)*z(j);
  lsum4(j E J): <<k E K| <<p E P| d(k,p)*y(j,k) >> >> =g= vl(j)*z(j);
  lsum5(k E K): <<j E J| y(j,k)>> =e= 1;
}}