function [ Jm,um,ym ] = plant(u,d,noise)


k1=d(1);
k2=d(2);
cain=2;
cbin=1.5;
V=500;
R=[1 0;0 1];


um=u;

x0=[0.0578
    0.0551
    0.7825];
u=u+noise(4:5);

Fa=u(1);
Fb=u(2);

opts=optimset('Display','off','TolX',1e-16,'TolFun',1e-16);
x=fsolve(@equations,x0,opts);
ca=x(1);
cb=x(2);
cc=x(3);

yn=[ 0.0578
    0.0551
    0.7825];
%measured variables
ym=[ca;cb;cc]+noise(1:3);



%cost function
Jm=(Fa+Fb)*cc*(Fa+Fb)*cc/(Fa*cain)-0.5*[Fa Fb]*R*[Fa;Fb];



%model equations at steady sate
    function sys=equations(x)
        
        ca=x(1);
        cb=x(2);
        cc=x(3);
        
        sys(1)=Fa*cain-(Fa+Fb)*ca-k1*ca*cb*V;
        sys(2)=Fb*cbin-(Fa+Fb)*cb-k1*ca*cb*V-2*k2*cb^2*V;
        sys(3)=-(Fa+Fb)*cc+k1*ca*cb*V;
        
    end



end

