var xA{1..3,0..5} >= 0;
var xp{1..3,1..5} >= 0;
param r1{1..3,1..5};
param r2{1..3,1..5};
param kA{1..3,1..5};
param k1{1..3,1..5};
param kp{1..3,1..5};
param kz{1..3,1..5};
param kL1{1..3,1..5};
param kL2{1..3,1..5};
param h{1..3,1..5};
param t{1..3,1..3};
param b{1..3,1..5};
minimize f:
-1*sum {m in 1..5} sum {j in 1..3} (r1[j,m]*xA[j,m] - r2[j,m]*xA[j,m]^2 - kA[j,m]*xA[j,m] - k1[j,m]*xp[j,m] - 
kp[j,m]*xp[j,m] - kz[j,m]*(xp[j,m]-xA[j,m-1])^2 - kL1[j,m]*(xp[j,m]-xA[j,m])^2 - kL2[j,m]*sum{s in 1..m} (xp[j,s]-xA[j,s]));

subject to cons1 {j in 1..3, m in 1..5}:
h[j,m] - xA[j,m] >= 0;

subject to cons2 {i in 1..3, m in 1..5}:
b[i,m] - sum {j in 1..3} t[i,j]*xp[j,m] >= 0;

subject to cons3 {j in 1..3, m in 1..5}:
sum {s in 1..m} (xp[j,s]-xA[j,s]) >= 0;

subject to cons4{j in 1..3}:
xA[j,0] = 0;

data;
var xA:=
1	1	80
2	1	100
3	1	400
1	2	100
2	2	200
3	2	200
1	3	100
2	3	250
3	3	400
1	4	50
2	4	200
3	4	500
1	5	50
2	5	200
3	5	500;

var xp:=
1	1	100
2	1	120
3	1	410
1	2	120
2	2	250
3	2	250
1	3	150
2	3	300
3	3	410
1	4	600
2	4	250
3	4	510
1	5	100
2	5	250
3	5	510;

param r1(tr):
	1		2		3:=
1	1000	520		910
2	1000	520		910
3	1000	520		1000
4	1100	600		1000
5	1100	600		1000;

param r2(tr):
	1	2	3:=
1	0.3	0.1	0.2	
2	0.3	0.1	0.2	
3	0.3	0.1	0.2	
4	0.3	0.1	0.2	
5	0.3	0.1	0.2;	

param kA(tr):
	1		2		3:=
1	120		65		105
2	150		65		105
3	150		80		120
4	170		80		120
5	170		80		120;

param k1(tr):
	1		2		3:=
1	150		75		140
2	150		75		140
3	150		75		140
4	170		90		150
5	170		90		150;

param kp(tr):
	1		2		3:=
1	160		75		140
2	160		75		140
3	160		75		140
4	180		90		150
5	180		90		150;

param kz(tr):
	1		2		3:=
1	0.02	0.01	0.015
2	0.2		0.1		0.15
3	0.25	0.1		0.15
4	0.25	0.15	0.15
5	0.25	0.15	0.15;

param kL1(tr):
	1		2		3:=
1	0.005	0.005	0.005
2	0.05	0.05	0.05
3	0.06	0.06	0.06
4	0.06	0.06	0.06
5	0.06	0.06	0.06;

param kL2(tr):
	1		2		3:=
1	80		45		75
2	80		45		75
3	100		45		90
4	100		50		90
5	100		50		90;

param h(tr):
	1		2		3:=
1	100		280		520
2	180		400		400
3	220		450		500
4	150		450		630
5	100		400		600;

param t(tr):
	1		2		3:=
1	0.6		0.4		0.1
2	0.3		0.1		0.12
3	0.36	0.08	0.06;

param b(tr):
	1		2		3:=
1	170		170		180
2	170		170		180
3	170		170		180
4	170		170		180
5	170		170		180;
solve;

display f;