********************************************************************
********************************************************************

This chapter is from the book:

Castillo, E., Conejo A.J., Pedregal, P., Garc\'{\i}a, R. and Alguacil, N. (2002). Building and Solving Mathematical Programming Models in Engineering and Science, Pure and Applied Mathematics Series, Wiley, New York.

Copyright © 2002 by John Wiley and Sons, Inc.

This material is used by permission of John Wiley and Sons, Inc.

********************************************************************
********************************************************************

$title Neural networks (neural2)

SETS
S number of data vectors/1*30/
I dimension of data input vectors/1*2/
J number of different learning functions/1*5/
R number of basic functions/1*3/;

ALIAS(S,S1)

PARAMETERS
X(I,S) input data vectors
Y(S) output data
powers(R) exponents of x in the basic functions
/1 0.5
2  1
3  2/;

X(I,S)=uniform(0,1)*0.5;
Y(S)=sqr(0.3+0.3*X('1',S)+0.7*X('2',S))+(uniform(0,1)-0.5)*0.01;

PARAMETERS
aux auxiliary parameter
x1 auxiliary parameter
x2 auxiliary parameter
x3 auxiliary parameter
f1 auxiliary parameter
f2 auxiliary parameter
f3 auxiliary parameter
ymin minimum value of Y(S)
ymax maximum value of Y(S)
acterror maximum prediction error
maxerror maximum allowed error for the bisection method;

maxerror=0.00001;

VARIABLES
z1 function1 to be optimized
z2 function2 to be optimized
z3 function3 to be optimized
z4 function21 to be optimized
z5 function21 to be optimized
W0 threshold value
W(I) weight associated with component i of input
alfa(R) coefficients of the neural function;

POSITIVE VARIABLES
epsilon error;

EQUATIONS
Q1 definition of z1
Q2 definition of z2
Q3 definition of z3
Q4 definition of z4
Q5 definition of z5
const1(S) upper error bound
const2(S) lower error bound
const3(S) upper error bound
const4(S) lower error bound
normalized normalized values
increasing(S,S1) the f function must be increasing;

Q1..z1=e=sum(S,sqr(Y(S)-arctan(W0+sum(I,W(I)*X(I,S)))));
Q2..z2=e=sum(S,sqr(W0+sum(I,W(I)*X(I,S))-sin(Y(S))/cos(Y(S))));
Q3..z3=e=epsilon;
Q4..z4=e=sum(S,sqr(W0+sum(I,W(I)*X(I,S))
            -sum(R,alfa(R)*Y(S)**(powers(R)))));
Q5..z5=e=epsilon;
const1(S)..W0+sum(I,W(I)*X(I,S))-epsilon=l=sin(Y(S))/cos(Y(S));
const2(S)..-W0-sum(I,W(I)*X(I,S))-epsilon=l=-sin(Y(S))/cos(Y(S));
const3(S)..W0+sum(I,W(I)*X(I,S))
      -sum(R,alfa(R)*Y(S)**(powers(R)))-epsilon=l=0.0;
const4(S)..-W0-sum(I,W(I)*X(I,S))
      +sum(R,alfa(R)*Y(S)**(powers(R)))-epsilon=l=0.0;
normalized..sum(R,alfa(R)*Y('1')**(powers(R)))=e=1;
increasing(S,S1)$(Y(S)<Y(S1))..sum(R,alfa(R)*Y(S)**(powers(R)))
      =l=sum(R,alfa(R)*Y(S1)**(powers(R)));;

MODEL reg1/Q1/;
MODEL reg2/Q2/;
MODEL reg3/Q3,const1,const2/;
MODEL reg4/Q4,normalized,increasing/;
MODEL reg5/Q5,const3,const4,normalized,increasing/;

file out/neural3.out/;
put out;

put "Data"/;
ymin=0.0;
ymax=0.0;
loop(S,
  loop(I,put X(I,S):12:3," & ";);
  put Y(S):12:3,"\\"/;
  if(ymin<Y(S),ymin=Y(S));
  if(ymax>Y(S),ymax=Y(S));
);
loop(J,
  if(ord(J) eq 1, SOLVE reg1 USING nlp MINIMIZING z1;
        put "z1=",z1.l:12:9/;);
  if(ord(J) eq 2, SOLVE reg2 USING nlp MINIMIZING z2;
        put "z2=",z2.l:12:9/;);
  if(ord(J) eq 3, SOLVE reg3 USING  lp MINIMIZING z3;
        put "z3=",z3.l:12:9/;);
  if(ord(J) eq 4, SOLVE reg4 USING nlp MINIMIZING z4;
        put "z4=",z4.l:12:9/;);
  if(ord(J) eq 5, SOLVE reg5 USING  lp MINIMIZING z5;
        put "z5=",z5.l:12:9/;);
  put "Weights:"/;
  put W0.l:12:3;
  loop(I,put " & ",W.l(I):12:3;);
  put " "/;
  if(((ord(J) eq 4) or (ord(J) eq 5)),
    put "Alfas:"/;
    loop(R,put " & ",alfa.l(R):12:3;);
    put " "/;
  );
  acterror=-10000;
  put "Data and fitted values:"//;

  loop(S,
    aux=W0.l+sum(I,W.l(I)*X(I,S));
    if(ord(J) le 3,
      f1=arctan(aux);
    else
      x1=ymin-(ymax-ymin)*0.1;
      x2=ymax+(ymax-ymin)*0.1;
      f1=sum(R,alfa.l(R)*abs(x1)**(powers(R)))-aux;
      f2=sum(R,alfa.l(R)*abs(x2)**(powers(R)))-aux;
      if(f1*f2>0.0,
        put "ERROR IN BISECTION "," S=",S.tl:3," aux=",
                aux:12:3," f1=",f1:12:3," f2=",f2:12:3/;
      else
        while(abs(x1-x2)>maxerror,
          x3=(x1+x2)*0.5;
          f3=sum(R,alfa.l(R)*x3**(powers(R)))-aux;
          if(f3*f1>0.0,
            x1=x3;f1=f3;
          else
            x2=x3;f2=f3;
          );
        );
      );
      f1=((x1+x2)*0.5);
    );
    loop(I,put X(I,S):12:3," & ";);
    aux=(f1-Y(S));
    if(abs(aux)>acterror,acterror=abs(aux));
    put Y(S):12:3," & ",aux:12:3"\\"/;
  );
  put "Maximum error=",acterror:12:9/;
);

********************************************************************
********************************************************************

$title Mesh

SETS
I number of elements/1*242/
J number of nodes per element/1*3/
K number of nodes/1*144/
N number of nodes per row/1*12/
D space dimension/1*2/
XFIXED(K) nodes with fixed X
YFIXED(K) nodes with fixed Y;

ALIAS(K,K1,K2);
\end{verbatim}
}\label{dynamicex1}{\footnotesize
\begin{verbatim}
PARAMETER
t number of squares per dimension
m number of elements per column
xmin minimum value of coordinate x
xmax maximum value of ccordinate x
ymin minimum value of coordinate y
ymax maximum value of ccordinate y
COORD(K,D) coordinates of the initial mesh nodes;

t=card(N);
xmin=-2;
xmax=2;
ymin=-2;
ymax=2;

m=2*t-2;

XFIXED(K)=no;
loop(N,
	XFIXED(K)$(ord(K) eq ord(N))=yes;
	XFIXED(K)$(ord(K) eq ord(N)+t*t-t)=yes;
	YFIXED(K)$(ord(K) eq t*ord(N))=yes;
	YFIXED(K)$(ord(K) eq (ord(N)-1)*t+1)=yes;
);

COORD(K,'1')=xmin+floor((ord(K)-1)/t)*(xmax-xmin)/(t-1);
COORD(K,'2')=ymin+mod(ord(K)-1,t)*(ymax-ymin)/(t-1);

PARAMETER

EL(I,J);

loop(I,
	if(mod(ord(I),2) eq 0,
		EL(I,'1')=(floor((ord(I)-1)/m))*t+floor(mod((ord(I)-1),m)/2+1);
		EL(I,'2')=(floor((ord(I)-1)/m))*t+floor(mod((ord(I)-1),m)/2+1)+1;
		EL(I,'3')=(floor((ord(I)-1)/m))*t+floor(mod((ord(I)-1),m)/2+1)+t+1;
	else
	 EL(I,'1')=(floor((ord(I)-1)/m))*t+floor(mod((ord(I)-1),m)/2+1);
		EL(I,'2')=(floor((ord(I)-1)/m))*t+floor(mod((ord(I)-1),m)/2+1)+t;
		EL(I,'3')=(floor((ord(I)-1)/m))*t+floor(mod((ord(I)-1),m)/2+1)+t+1;
	);
);

VARIABLES
z1 function to be optimized
X(K) coordinate of node K
Y(K) coordinate of node K
Z(K) coordinate of node K
area(I) area of element I
mean mean area of elements;

X.l(K)=COORD(K,'1');
Y.l(K)=COORD(K,'2');

X.fx(XFIXED)=COORD(XFIXED,'1');
Y.fx(YFIXED)=COORD(YFIXED,'2');

EQUATIONS
Q objective function to be optimized
Carea(I)
Cf(K) points are in surface
Cmean;

Q..z1=e=sum(I,sqr(area(I)-mean))/card(I);

Carea(I)..area(I)=e=sum((K,K1,K2)$((EL(I,'1') 
   =ord(K)) and (EL(I,'2')=ord(K1)) and (EL(I,'3')=ord(K2))),
sqr((Y(K1)-Y(K))*(Z(K2)-Z(K))-(Y(K2)-Y(K))*(Z(K1)-Z(K)))+
sqr((Z(K1)-Z(K))*(X(K2)-X(K))-(Z(K2)-Z(K))*(X(K1)-X(K)))+
sqr((X(K1)-X(K))*(Y(K2)-Y(K))-(X(K2)-X(K))*(Y(K1)-Y(K))));
Cf(K)..Z(K)=e=(X(K)*Y(K)*(sqr(X(K))-sqr(Y(K))))/(sqr(X(K))+sqr(Y(K)));
Cmean..mean=e=sum(I,area(I))/card(I);

MODEL mallas/ALL/;

SOLVE mallas USING nlp MINIMIZING z1;

file out/mallas3.out/;
put out;

if((card(I) ne (2*(t-1)*(t-1))),
	put "Error, card(I) must be equal to ",(2*(t-1)*(t-1)):8:0/;
else
	if((card(K) ne t*t),
		put "Error, card(K) must be equal to ",(t*t):8:0/;
	else
put "z1=",z1.l:12:8," mean=",mean.l:12:8, " cpu=",mallas.resusd:9:2,
     " modelstat=",mallas.modelstat:3:0," solvestat=",mallas.solvestat:3:0/;

loop(K,
	put "X(",K.tl:3,")=",COORD(K,'1'):12:8," Y(",K.tl:3,")=",COORD(K,'2'):12:8/;
);

loop(I,
	put "I=",i.tl:3," Nodes=",EL(I,'1'):5:0," ",EL(I,'2'):5:0," ",EL(I,'3'):5:0/;
);

loop(K,
	put "K=",K.tl:3," X=",X.l(K):12:6," Y=",Y.l(K):12:6," Z=",Z.l(K):12:6/;
);

loop(I,
	put "I=",I.tl:3," area=",area.l(I):12:6/;
);

put "COORDINATES"/;
put "{";

loop(K,
	if(ord(K)<card(K),
		put "{",X.l(K):12:8,",",Y.l(K):12:8,",",Z.l(K):12:8,"},"/;
	else
		put "{",X.l(K):12:8,",",Y.l(K):12:8,",",Z.l(K):12:8,"}}"/;
	);
);

put "ELEMENTS"/;
put "{";

loop(I,
	if(ord(I)<card(I),
		put "{",EL(I,'1'):5:0,",",EL(I,'2'):5:0,",",EL(I,'3'):5:0,"},"/;
	else
		put "{",EL(I,'1'):5:0,",",EL(I,'2'):5:0,",",EL(I,'3'):5:0,"}}"/;
	););););
\end{verbatim}

********************************************************************
********************************************************************

$title Compatibility Method I

SETS
I number of rows/1*6/
J number of columns/1*6/;

ALIAS(I,I1);
ALIAS(J,J1);

PARAMETER
PP(I,J) Joint probability matrix used to generate conditionals
A(I,J) conditional probability of X given Y
B(I,J) conditional probability of Y given X;

PP(I,J)=round(uniform(0,5));
PP(I,J)=PP(I,J)/sum((I1,J1),PP(I1,J1));
A(I,J)=PP(I,J)/sum(I1,PP(I1,J));
B(I,J)=PP(I,J)/sum(J1,PP(I,J1));

VARIABLE
z value of the objective function;

POSITIVE VARIABLES
P(I,J) Joint probability looked for
TAU(I) Marginal probability of X
ETA(J) Marginal probability of Y;

EQUATIONS
ZDEF1 Function to be optimized in model 1
ZDEF2 Function to be optimized in model 2
ZDEF3 Function to be optimized in model 3
CONST11(I,J)
CONST12(I,J)
CONST2(I,J)
CONST3(I,J)
NORM1 P add up to one
NORM2 TAU add up to one
NORM3 ETA add up to one;

ZDEF1..z=e=P('1','1');
ZDEF2..z=e=TAU('1');
CONST11(I,J)..P(I,J)-A(I,J)*sum(I1,P(I1,J))=e=0;
CONST12(I,J)..P(I,J)-B(I,J)*sum(J1,P(I,J1))=e=0;
CONST2(I,J)..ETA(J)*A(I,J)-TAU(I)*B(I,J)=e=0;
CONST3(I,J)..A(I,J)*sum(I1,TAU(I1)*B(I1,J))-TAU(I)*B(I,J)=e=0;
NORM1..sum((I,J),P(I,J))=e=1;
NORM2..sum(I,TAU(I))=e=1;
NORM3..sum(J,ETA(J))=e=1;

MODEL model1/ZDEF1,CONST11,CONST12,NORM1/;
MODEL model2/ZDEF2,CONST2,NORM2,NORM3/;
MODEL model3/ZDEF2,CONST3,NORM2/;

file out/comp1.out/;
put out;

put "Initial Matrix PP:"/;
loop(I,loop(J,put PP(I,J):7:4;);put ""/;);
put "Matrix A:"/;
loop(I,loop(J,put A(I,J):7:4;);put ""/;);
put "Matrix B:"/;
loop(I,loop(J,put B(I,J):7:4;);put ""/;);

SOLVE model1 USING lp MINIMIZING z;
put "Final Matrix P Obtained by Model 1:"/; 
loop(I,loop(J,put P.l(I,J):7:4;);put ""/;);
put "z=",z.l," modelstat=",model1.modelstat," solvestat=",model1.solvestat/;
SOLVE model2 USING lp MINIMIZING z;
put "Final Matrix P Obtained by Model 2:"/;
loop(I,loop(J,put (TAU.l(I)*B(I,J)):7:4;);put ""/;);
put "z=",z.l," modelstat=",model2.modelstat," solvestat=",model2.solvestat/;
SOLVE model3 USING lp MINIMIZING z;
put "Final Matrix P Obtained by Model 3:"/;
loop(I,loop(J,put (TAU.l(I)*B(I,J)):7:4;);put ""/;);
put "z=",z.l," modelstat=",model3.modelstat," solvestat=",model3.solvestat/;

********************************************************************
********************************************************************

$title Compatibility Methods 1, 2 and 3

SETS
I number of rows/1*4/
J number of columns/1*4/;

ALIAS(I,I1);
ALIAS(J,J1);

PARAMETER
PP(I,J) Joint probability matrix used to generate conditionals
A(I,J) conditional probability of X given Y
B(I,J) conditional probability of Y given X;

PP(I,J)=round(uniform(0,5));
PP(I,J)=PP(I,J)/sum((I1,J1),PP(I1,J1));
A(I,J)=PP(I,J)/sum(I1,PP(I1,J));
B(I,J)=PP(I,J)/sum(J1,PP(I,J1));

* The line below modifies the B matrix to get incompatibility

VARIABLE
z value of the objective function
epsilon Maximum error in joint probability;

POSITIVE VARIABLES
P(I,J) Joint probability looked for
TAU(I) Marginal probability of X
ETA(J) Marginal probability of Y;

EQUATIONS
ZDEF1 Function to be optimized in model 1
ZDEF2 Function to be optimized in model 2
ZDEF3 Function to be optimized in model 3
CONST111(I,J)
CONST112(I,J)
CONST121(I,J)
CONST122(I,J)
CONST21(I,J)
CONST22(I,J)
CONST31(I,J)
CONST32(I,J)
NORM1 P add up to one
NORM2 TAU add up to one
NORM3 ETA add up to one;

ZDEF1..z=e=epsilon;
CONST111(I,J)..P(I,J)-A(I,J)*sum(I1,P(I1,J))=l=epsilon;
CONST112(I,J)..-epsilon=l=P(I,J)-A(I,J)*sum(I1,P(I1,J));
CONST121(I,J)..P(I,J)-B(I,J)*sum(J1,P(I,J1))=l=epsilon;
CONST122(I,J)..-epsilon=l=P(I,J)-B(I,J)*sum(J1,P(I,J1));
CONST21(I,J)..ETA(J)*A(I,J)-TAU(I)*B(I,J)=l=epsilon;
CONST22(I,J)..-epsilon=l=ETA(J)*A(I,J)-TAU(I)*B(I,J);
CONST31(I,J)..A(I,J)*sum(I1,TAU(I1)*B(I1,J))-TAU(I)*B(I,J)=l=epsilon;
CONST32(I,J)..-epsilon=l=A(I,J)*sum(I1,TAU(I1)*B(I1,J))-TAU(I)*B(I,J);
NORM1..sum((I,J),P(I,J))=e=1;
NORM2..sum(I,TAU(I))=e=1;
NORM3..sum(J,ETA(J))=e=1;

MODEL model1/ZDEF1,CONST111,CONST112,CONST121,CONST122,NORM1/;
MODEL model2/ZDEF1,CONST21,CONST22,NORM2,NORM3/;
MODEL model3/ZDEF1,CONST31,CONST32,NORM2/;

file out/comp1.out/;
put out;

put "Initial Matrix PP:"/;
loop(I,loop(J,put PP(I,J):7:4," & ";);put ""/;);
put "Matrix A:"/;
loop(I,loop(J,put A(I,J):7:4," & ";);put ""/;);
put "Matrix B:"/;
loop(I,loop(J,put B(I,J):7:4," & ";);put ""/;);

SOLVE model1 USING lp MINIMIZING z;
put "Final Matrix P Obtained by Model 1:"/; 
loop(I,loop(J,put P.l(I,J):7:4," & ";);put ""/;);
put "z=",z.l:12:8," modelstat=",model1.modelstat," solvestat=",model1.solvestat/;
SOLVE model2 USING lp MINIMIZING z;
put "Final Matrix P Obtained by Model 2:"/;
loop(I,loop(J,put (TAU.l(I)*B(I,J)):7:4," & ";);put ""/;);
put "z=",z.l:12:8," modelstat=",model2.modelstat," solvestat=",model2.solvestat/;
SOLVE model3 USING lp MINIMIZING z;
put "Final Matrix P Obtained by Model 3:"/;
loop(I,loop(J,put (TAU.l(I)*B(I,J)):7:4," & ";);put ""/;);
put "z=",z.l:12:8," modelstat=",model3.modelstat," solvestat=",model3.solvestat/;

********************************************************************
********************************************************************

$ title Regression

SETS
I number of data points /1*30/
P Number of parameters/1*5/
J number of regression models/1*3/
SS(I) subset of data points used in analysis;

PARAMETER
C(P) regression coefficients used in the simulation
/1 0.2
2 0.1
3 -0.3
4 0.4
5 -0.2/;

PARAMETER
Z(I,P) observed and predicted variables;

Z(I,P)=uniform(0,1);
Z(I,'1')=C('1')+sum(P$(ord(P)>1),C(P)*Z(I,P))+uniform(0,1)*0.05;
Z('7','1')=C('1')+sum(P$(ord(P)>1),C(P)*Z('7',P))-0.05*2;
Z('15','1')=C('1')+sum(P$(ord(P)>1),C(P)*Z('15',P))+0.05*2;
Z('16','1')=C('1')+sum(P$(ord(P)>1),C(P)*Z('16',P))-0.05*2;

PARAMETER
Y(I) Observed values
Y1(I) Fitted values with all data points
Y2(I) Fitted values with one data point removed
z0 auxiliary value;

Y(I)=Z(I,'1');

POSITIVE VARIABLES
EPSILON1(I) error associated with data I
EPSILON error;

VARIABLES
z1 objective function value
z2 objective function value
z3 objective function value
BETA(P);

EQUATIONS
dz1 objective function 1 value definition
dz2 objective function 2 value definition
dz3 objective function 3 value definition
lower1(I) lower bound of error
upper1(I) upper bound of error
lower2(I) lower bound of error
upper2(I) upper bound of error;

dz1..z1=e=sum(SS,EPSILON1(SS));
dz2..z2=e=EPSILON;
dz3..z3=e=sum(SS,sqr(Y(SS)-sum(P,Z(SS,P)*BETA(P))));
lower1(SS)..Y(SS)-sum(P,Z(SS,P)*BETA(P))=l=EPSILON1(SS);
upper1(SS)..-Y(SS)+sum(P,Z(SS,P)*BETA(P))=l=EPSILON1(SS);
lower2(SS)..Y(SS)-sum(P,Z(SS,P)*BETA(P))=l=EPSILON;
upper2(SS)..-Y(SS)+sum(P,Z(SS,P)*BETA(P))=l=EPSILON;

MODEL regreslav/dz1,lower1,upper1/;
MODEL regresMM/dz2,lower2,upper2/;
MODEL regresLS/dz3/;

file out/regress1.out/;
put out;

put "Data"/;

loop(I,
	put I.tl:3;
	loop(P,put " & ",Z(I,P):8:4;);
	put "\\"/;
);
Z(I,'1')=1.0;

loop(J,
  SS(I)=yes;
  if(ord(J) eq 1,SOLVE regreslav USING lp MINIMIZING z1;Z0=z1.l);
  if(ord(J) eq 2,SOLVE regresMM USING lp MINIMIZING z2;Z0=z2.l);
  if(ord(J) eq 3,SOLVE regresLS USING nlp MINIMIZING z3;Z0=z3.l);
  put "All data points z0=",z0:12:8 ;
    loop(P,put " & ",BETA.l(P):10:5;);
    put "\\"/;
    Y1(I)=sum(P,Z(I,P)*BETA.l(P));
  loop(I,
    SS(I)=no;
    if(ord(J) eq 1,SOLVE regreslav USING lp MINIMIZING z1;Z0=z1.l);
    if(ord(J) eq 2,SOLVE regresMM USING lp MINIMIZING z2;Z0=z2.l);
    if(ord(J) eq 3,SOLVE regresLS USING nlp MINIMIZING z3;Z0=z3.l);
    put "Removing point ",I.tl:3," z0=",z0:12:8;
    loop(P,put " & ",BETA.l(P):10:5;);
    put "\\"/;
    Y2(I)=sum(P,Z(I,P)*BETA.l(P))-Y1(I);
    SS(I)=yes;
  );
  loop(I,put I.tl:3,Y2(I):10:6/;);
);

********************************************************************
********************************************************************

$title FUNCTIONAL1 n=20
SET J /0*20/;
VARIABLES z,x(J);
  x.fx(J)$(ord(J) eq 1) = 0;
  x.fx(J)$(ord(J) eq card(J)) = 1;
SCALAR n;
  n = card(J)-2;
EQUATION
  cost   objective function;
cost.. z =e= SUM(J$(ord(J) lt card(J)),
(sqrt(ord(J))-sqrt(ord(J)-1))*sqrt(1+sqr(n+1)*sqr(x(J+1)-x(J))));
MODEL funct1 /all/;
SOLVE funct1 USING nlp MINIMIZING z;

********************************************************************
********************************************************************

$title FUNCTIONAL2 n=10
SET J /0*10/;
VARIABLES z,x(J);
  x.fx(J)$(ord(J) eq 1) = 0;
  x.fx(J)$(ord(J) eq card(J)) = 0;
  x.l(J)=-0.5;
SCALAR n;
  n = card(J)-2;
EQUATION
  cost   objective function
  rest   equality constraint;
cost.. z =e= SUM(J$(ord(J) lt card(J)),
                  sqrt(1+sqr(n+1)*sqr(x(J+1)-x(J)))*(x(J+1)+x(J)));
rest.. (n+1)*1.3116 =e= SUM(J$(ord(J) lt card(J)),sqrt(1+sqr(n+1)*sqr(x(J+1)-x(J))));
MODEL funct2 /all/;
SOLVE funct2 USING nlp MINIMIZING z;

********************************************************************
********************************************************************

$title FUNCTIONAL3 n=30
SET K /1*30/;
VARIABLES z,u(K);
  u.up(K) = 1;
  u.lo(K) = 0;
SCALAR n;
  n = card(K);
EQUATION
  cost   objective function
  rest   equality constraint;
cost.. z =e= SUM(K,sqr(u(K)));
rest.. 3+(5/6) =e= (9/(2*sqr(n)))*SUM(K,(2*n-2*ord(K)+1)*u(K));
MODEL funct3 /all/;
SOLVE funct3 USING nlp MINIMIZING z;

********************************************************************
********************************************************************

$title FUNCTIONAL n=20, azero
SET J /1*20/;
SCALAR pi the pi number /3.1416/;
VARIABLES z,theta,u(J);
 u.lo(J)=-1;
 u.up(J)=1;
 theta.lo=0;
 theta.up=pi/2;
 theta.l=pi/4;
 u.l(J)=0.;
SCALARS n, azero;
  n = card(J);
  alpha=-7;
EQUATION
  cost   objective function
  const1
  const2;
cost.. z =e= SIN(theta)/COS(theta);;
const1..
alpha*COS(n*theta)*POWER((1+POWER((SIN(theta)/COS(theta)),2)),n)+(SIN(theta)
/COS(theta))*
         SUM(J,u(J)*SIN((n-ORD(J))*theta)*
           POWER((1+POWER((SIN(theta)/COS(theta)),2)),n-ORD(J))
         ) =E= 0;
const2..
-azero*SIN(n*theta)*POWER((1+POWER((SIN(theta)/COS(theta)),2)),n)+(SIN(theta
)/COS(theta))*
         SUM(J,u(J)*COS((n-ORD(J))*theta)*
         POWER((1+POWER((SIN(theta)/COS(theta)),2)),n-ORD(J))
         ) =E= 0;
MODEL funct5 /all/;
SOLVE funct5 USING nlp MINIMIZING z;

********************************************************************
********************************************************************

$title TRAFFIC ASSIGNMENT PROBLEM.

** Sets are declared in first place.
** Set N is the set of nodes of the road network.
** Set A(N,N) is the set of links.
** Set W is the set of O--D pairs

SET
N set of nodes /I1*I13/
W pairs        /W1*W4/
A(N,N) set of links
/I1.(I5,I12)
I4.(I5,I9)
I5.(I6,I9)
I6.(I7,I10)
I7.(I8,I11)
I8.(I2)
I9.(I10,I13)
I10.(I11)
I11.(I2,I3)
I12.(I6,I8)
I13.(I3)/

ODE(N,W) set of origins by demand
/I1.(W1,W2)
I4.(W3,W4)/

DDE(N,W) set of destinations by demand
/I2.(W1,W3)
I3.(W2,W4)/;

** The set of nodes I must be  duplicated to refer to its different element
*  in the same constraint.

ALIAS(N,I)
ALIAS(N,J)

PARAMETER
G(W) demand in the O--D pair  W
/W1 400
W2 800
W3 600
W4 200/;

TABLE  CDATA(N,N,*)  link cost parameters
          C0      D       M
 I1.I5    7.0  0.00625   2.0
 I1.I12   9.0  0.00500   2.0
 I4.I5    9.0  0.00500   2.0
 I4.I9   12.0  0.00250   2.0
 I5.I6    3.0  0.00375   2.0
 I5.I9    9.0  0.00375   2.0
 I6.I7    5.0  0.00625   2.0
 I6.I10  13.0  0.00250   2.0
 I7.I8    5.0  0.00625   2.0
 I7.I11   9.0  0.00625   2.0
 I8.I2    9.0  0.00625   2.0
 I9.I10  10.0  0.00250   2.0
 I9.I13   9.0  0.00250   2.0
 I10.I11  6.0  0.00125   2.0
 I11.I2   9.0  0.00250   2.0
 I11.I3   8.0  0.00500   2.0
 I12.I6   7.0  0.00125   2.0
 I12.I8  14.0  0.00500   2.0
 I13.I3  11.0  0.00500   2.0;

** Optimization variables are declared.

VARIABLES
z         objective function variable
f(I,J)    the flow at the link I-J
fc(W,I,J) the flow of the commodity W at the link I-J
cos(I,J)  cost in the link I-J at equilibrium;

POSITIVE VARIABLE fc(W,N,N);

** Constraints are declared.

EQUATIONS
COSTZ Objective function
BALANCE(W,I) conservation of flow condition for the commodity W.
FLOW(I,J) total link flow
COST(I,J) cost of the link I-J;

** The objective function is the sum of the integrated cost
** in all links.
** The equation FLOW illustrate the use of dollar operation
** to restrict the number of constraints generated to less
** than implied by the domain of the defining sets.

COSTZ.. z=E= SUM((I,J)$A(I,J),CDATA(I,J,'C0')*f(I,J)
              +CDATA(I,J,'D')*f(I,J)**CDATA(I,J,'M'));
BALANCE(W,I) .. SUM(J$A(I,J),fc(W,I,J))-SUM(J$A(J,I),fc(W,J,I))
              =E=G(W)$ODE(I,W)-G(W)$DDE(I,W);
FLOW(I,J)$A(I,J)..SUM(W,fc(W,I,J))=E=f(I,J);
COST(I,J)$A(I,J)..cos(I,J)=E=CDATA(I,J,'C0')
             +2*CDATA(I,J,'D')*f(I,J);

** The next two sentences define the ND (Nguyen-Dupis) problem,
** considering all the above constraints, and direct GAMS to
** solve the problem using the nlp solver.

MODEL ND /ALL/;
SOLVE ND USING nlp MINIMIZING z;

********************************************************************
********************************************************************

$title TRAFFIC ASSIGNMENT PROBLEM WITH ELASTIC DEMAND.

** Sets are declared in first place.
** Set N is the set of nodes of the road network.
** Set A is the set of links.
** Set W is the set of O--D pairs

SET
N set of nodes /I1*I13/
W pairs         /W1*W4/
A(N,N) set of links
/I1.(I5,I12)
I4.(I5,I9)
I5.(I6,I9)
I6.(I7,I10)
I7.(I8,I11)
I8.(I2)
I9.(I10,I13)
I10.(I11)
I11.(I2,I3)
I12.(I6,I8)
I13.(I3)/

ODE(N,W) set of origins by demand
/I1.(W1,W2)
I4.(W3,W4/

DDE(N,W) set of destinations by demand
/I2.(W1,W3)
I3.(W2,W4)/;

ALIAS(N,I)
ALIAS(N,J)

PARAMETER
G(W) demand in O--D pair  W
/W1 600
W2 1000
W3 800
W4 400/

Cb(W) travel cost by public transport in O--D pair  W
/W1 41
W2 35
W3 43
W4 40/

BETA
ALPHAa
ALPHAb;
BETA=0.1;
ALPHAa=-2.;
ALPHAb=0.;

TABLE  CDATA(N,N,*)  link cost parameters
          C0      D
 I1.I5    7.0  0.00625
 I1.I12   9.0  0.005
 I4.I5    9.0  0.005
 I4.I9   12.0  0.0025
 I5.I6    3.0  0.00375
 I5.I9    9.0  0.00375
 I6.I7    5.0  0.00625
 I6.I10  13.0  0.0025
 I7.I8    5.0  0.00625
 I7.I11   9.0  0.00625
 I8.I2    9.0  0.00625
 I9.I10  10.0  0.0025
 I9.I13   9.0  0.0025
 I10.I11  6.0  0.00125
 I11.I2   9.0  0.0025
 I11.I3   8.0  0.005
 I12.I6   7.0  0.00125
 I12.I8  14.0  0.005
 I13.I3  11.0  0.005;


** Optimization variables are declared.

VARIABLES
 z         objective function variable
 f(N,N)    flow link
 fc(W,N,N) is the flow link of the commodity W
 cos(N,N)  cost in the link at equilibrium
 ga(W)     demand by car
 gb(W)     demand by public transport;

POSITIVE VARIABLE fc(W,N,N);

ga.LO(W)=0.01;
gb.LO(W)=0.01;

** Constraints are declared.

EQUATIONS
 COST Objective function
 BALANCE(W,I) conservation of flow condition for the commodity W.
 FLOW(I,J) total link flow
 MODAL(W) Modal split of the demand
 COSTE(I,J) cost at equilibrium;

COST .. z=E= SUM((I,J)$A(I,J),CDATA(I,J,'C0')*f(I,J)
             +CDATA(I,J,'D')*f(I,J)**CDATA(I,J,'M'))
+SUM(W,Cb(W)*gb(W))+1/BETA* SUM(W,ga(W)*(-1+ALPHAa+LOG(ga(W)))
             + gb(W)*(-1+ALPHAb+LOG(gb(W))) );
BALANCE(W,I) .. SUM(J$A(I,J),fc(W,I,J))-SUM(J$A(J,I),fc(W,J,I))
                =E=ga(W)$ODE(I,W)-ga(W)$DDE(I,W);
FLOW(I,J)$A(I,J)..SUM(W,fc(W,I,J))=E=f(I,J);
COSTE(I,J)$A(I,J)..cos(I,J)=E=CDATA(I,J,'C0')
                              +2*CDATA(I,J,'D')*f(I,J);
MODAL(W).. G(W)=E=ga(W)+gb(W);
MODEL nd /ALL/;
SOLVE nd USING nlp MINIMIZING z;

********************************************************************
********************************************************************

$title TRAFFIC DISTRIBUTION-ASSIGNMENT PROBLEM.

** Sets are declared in first place.
** Set N is the set of nodes of the road network.
** Set A is the set of links.
** Set W is the set of O--D pairs

SET
 N set of nodes /I1*I13/
 W pairs        /W1*W4/
 O origins      /I1
                 I4/
 D destinations /I2
                 I3/
OW(O,W)
/I1.(W1,W2)
I4.(W3,W4)/
DW(D,W)
/I2.(W1,W3)
I3.(W2,W4)/

A(N,N) set of links
/I1.(I5,I12)
I4.(I5,I9)
I5.(I6,I9)
I6.(I7,I10)
I7.(I8,I11)
I8.(I2)
I9.(I10,I13)
I10.(I11)
I11.(I2,I3)
I12.(I6,I8)
I13.(I3)/

ODE(N,W) set of origins by demand
/I1.W1
I1.W2
I4.W3
I4.W4/
DDE(N,W) set of destinations by demand
/I2.W1
I3.W2
I2.W3
I3.W4/;

ALIAS(N,I)
ALIAS(N,J)

PARAMETER
Oi(O) number of trips emanating from O
/I1 1200
I4 800/

Dj(D) number of trips attracted to D
/I2 1000
I3 1000/

BETA;
BETA=0.2;

TABLE  CDATA(N,N,*)  link cost parameters
          C0      D       M
 I1.I5    7.0  0.00625   2.0
 I1.I12   9.0  0.005     2.0
 I4.I5    9.0  0.005     2.0
 I4.I9   12.0  0.0025    2.0
 I5.I6    3.0  0.00375   2.0
 I5.I9    9.0  0.00375   2.0
 I6.I7    5.0  0.00625   2.0
 I6.I10  13.0  0.0025    2.0
 I7.I8    5.0  0.00625   2.0
 I7.I11   9.0  0.00625   2.0
 I8.I2    9.0  0.00625   2.0
 I9.I10  10.0  0.0025    2.0
 I9.I13   9.0  0.0025    2.0
 I10.I11  6.0  0.00125   2.0
 I11.I2   9.0  0.0025    2.0
 I11.I3   8.0  0.005     2.0
 I12.I6   7.0  0.00125   2.0
 I12.I8  14.0  0.005     2.0
 I13.I3  11.0  0.005     2.0;

** Optimization variables are declared.

VARIABLES
 z         objective function variable
 f(N,N)    flow link
 fc(W,N,N) is the flow link of the commodity W
 cos(N,N)  cost in the link at equilibrium
 g(W)      demand in the pair W;

POSITIVE VARIABLE fc(W,N,N);

g.LO(W)=0.01;

** Constraints are declared.

EQUATIONS
 COST Objective function
 BALANCE(W,I) conservation of flow condition for the commodity W.
 FLOW(I,J) total link flow
 ORIGIN(O)
 DESTIN(D)
 COSTE(I,J) cost at equilibrium;

COST .. z=E= SUM((I,J)$A(I,J),CDATA(I,J,'C0')*f(I,J)
             +CDATA(I,J,'D')*f(I,J)**CDATA(I,J,'M'))
             +(1/BETA)*SUM(W,g(W)*(LOG(g(W))-1 ));
BALANCE(W,I) .. SUM(J$A(I,J),fc(W,I,J))-SUM(J$A(J,I),fc(W,J,I))
                =E=g(W)$ODE(I,W)-g(W)$DDE(I,W);
FLOW(I,J)$A(I,J)..SUM(W,fc(W,I,J))=E=f(I,J);
COSTE(I,J)$A(I,J)..cos(I,J)=E=CDATA(I,J,'C0')
                              +2*CDATA(I,J,'D')*f(I,J);
ORIGIN(O)..Oi(O)=E= SUM(W,g(W)$OW(O,W));
DESTIN(D)..Dj(D)=E= SUM(W,g(W)$DW(D,W));
MODEL nd /ALL/;
SOLVE nd USING nlp MINIMIZING z;