function q = Q_P(Kp, G,t)
PID = tf([Kp],[1]);
C = feedback(G*PID,1);
[y,T] = step(C,t);
Es = zeros(rows(y),1);
for i=1:rows(y)
    Es(i)=(1-y(i))*(1-y(i));
end
q = trapz(T, Es);
endfunction

function qp = Q_PID(Kp)
G = tf([2], [2 3 1]);
t = 0:0.1:10;
qp = Q_P(Kp, G, t);
endfunction

fminsearch(@Q_PID,0)
ans =  364.75