function [x,u,K,S]=scaopt(a,b,q,r,s,x0,N)
S(N+1)=s;
for k=N:-1:1
    K(k)=(a*b*s)/(r+s*b^2);
    s=q+(r*s*a^2)/(r+s*b^2);
    S(k)=s;
end

x(1)=x0;
for k=1:N
    u(k)=-K(k)*x(k);
    x(k+1)=a*x(k)+b*u(k);
end