function [DeltaJ] = eval_gradient(A,B,K,Sigma_0,Q,R)
%Compute the gradient

P = dlyap((A-B*K)',Q+K'*R*K);
Sigma = dlyap(A-B*K,Sigma_0);
DeltaJ = 2*((R+B'*P*B)*K-B'*P*A)*Sigma;


end