function Ceq = col_eqconstr(u1,xb)

global d u NT

u = [u1 ; d(1:3)] ; %L, V, F, zF, qF

x0 = 0.5*ones(1,82)' ;
options = optimset('Display','none','TolFun',1e-10) ;
x = fsolve(@cola_lv,x0,options) ;

xB = x(1) ; 
xD = x(NT) ;

Ceq = [xB-xb ;
    0.95-xD] ;

