Help with modeling this problem

Hello everybody,
I have modeled the following problem. However, whenever I run my program with certain data (I am sure they are correct and appropriate), it does not give me what I expect. I guess there must be some problems with the way I code it. Any help will be appreciated.
model.JPG
This the code I have written. I even wrote R_q ^a and K_q ^a in different ways but I still have problems with defining these two sets. GAMS gives some errors concerning the variable “x”.

Sets
k /1*25/

q /1*86/

a /1*86/
;

alias(i,k);


Variable
o;

Binary Variables
x(k)
y(a,q,i)
z(q);

Parameters
p
f_ev(q)
d(q)
f_phev(q)
d(q);

p = 3;

call gdxxrw Data.xlsx par d rng=edges!A1:B86 rdim=1 cdim=0 gdxin Data.gdx
load d gdxin

call gdxxrw Data.xlsx par f_ev rng=flows!A1:B86 rdim=1 cdim=0 gdxin Data.gdx
load f_ev gdxin

call gdxxrw Data.xlsx par f_phev rng=flows!A1:B86 rdim=1 cdim=0 gdxin Data.gdx
load f_phev gdxin


Equations
Obj
Const1(a,q,i)
Const2(a,q)
Const3(q)
Const4;

Obj…
o =e= sum(q,f_ev(q)*d(q)*z(q)) + sum((a,q,i),f_phev(q)*d(q)*y(a,q,i));

Const1(a,q,i)…
y(a,q,i) =l= x(i);

Const2(a,q)…
sum(i,y(a,q,i)) =l= 1;

Const3(q)…
z(q) =l= sum(k,x(k));

Const4…
sum(k,x(k)) =e= p;

model CSLM_PHEV /all/;
solve CSLM_PHEV maximizing o using MIP;
display x.l,y.l,z.l,o.l,f_ev,f_phev,d;

Hi Parisa
If you would like to raise the probability of an answer on your question from 0 to 99%, you should follow the guidelines to ask your question (or first do a search).
Cheers
Renger