function J = twocol_dyn_obj(u)

global d price x_steady config

%Calculate objective function for two columns by running a dynamic model to
%steady state with composition control in top and bottom. Inputs u are the
%set points for compositions.

%Simulate for 5000 seconds

[tsol,xsol] = ode15s(@twocol_dyn,[0 5000],x_steady,[],u,d);

%Take x vector at final time step and retrieve necessary values from model:

%[xprime D1 B1 L1 V1 xD1 xB1 e1 D2 B2 L2 V2 xD2 xB2 e2] = twocol_dyn(xsol(end,:),u,d) ;
[xprime D1 B1 L1 V1 xD1 xB1 e1 D2 B2 L2 V2] = twocol_dyn(xsol(end,:),u,d);  %Don't need xD2, xB2 and e2 for now.
if max(xprime) > 1e-6
    [tsol,xsol] = ode15s(@twocol_dyn,[5000 25000],xsol(end,:),[],u,d);
    [xprime D1 B1 L1 V1 xD1 xB1 e1 D2 B2 L2 V2] = twocol_dyn(xsol(end,:),u,d);  %Don't need xD2, xB2 and e2 for now.
    if max(xprime) > 1e-6
        disp('Steady state not reached after 25000 seconds')
    end
end

%Now calculate objective function:
pF=price(1) ; pV = price(2) ; pA = price(3) ; pB = price(4) ; pC = price(5) ;

J = d(1)*pF + pV*(V1 + V2) - D1*pA - D2*pB - B2*pC ;