Observer-based feedback for undamped, unstable oscillator; simulate combined system
Problem 4.9

• revised, Oct. 24, 2022:
  – simplify use of EstimatorRegulator command for closed loop
  – miscellaneous code clean-up

Clear["Global`*"];
G0=(1/((s^2)-1)); G0tf=TransferFunctionModel[G0,s];
G0ss = StateSpaceModel[G0tf]; 
K = StateFeedbackGains[G0ss,{-1,-1}]   (* Feedback gains *);
L=EstimatorGains[G0ss,{-2,-2}]         (* Observer gains; note sign convention! *);

{A,b,c,d}=Normal[G0ss]  (* get A, b, c matrices *);
A1=A-b.K;
A2 = A-b.K-L.c; Ac = ArrayFlatten[(({{A, -b.K}, {L.c, A2}}))]; (* for "big" system *)
kr = ((-1)/(c.Inverse[A1].b));             (* feedforward gain for reference command *)
Br= ArrayFlatten[(({{b.kr}, {b.kr}}))];  (* input coupling for reference [not used here] *)
Bd =  ArrayFlatten[(({{b}, {0}, {0}}))];  (* assume disturbance couples through input *)
Btot = ArrayFlatten[ (({{Br, Bd}}))];   (* assemble the "big" input-coupling matrix *)
G1ss = StateSpaceModel[{Ac,Btot}]; (* the full system, w/o output coupling *)
u=(({{0}, {DiracDelta[t]}}))  (* no reference; disturb state *);

{θ,θd,θh,θhd} = StateResponse[G1ss,u,t]//Simplify;
tmax=10; 
p1=Plot[{θ,θh},{t,0,tmax},PlotRange->Full,PlotStyle ->{Thick},ImageSize->250];
p2=Plot[{θd,θhd},{t,0,tmax},PlotRange->Full,PlotStyle ->{Thick},ImageSize->250];
Grid[{{p1,p2}},Spacings->{4}]

Or do all this directly using a special Mathematica system command:

closedLoop=EstimatorRegulator[G0ss,{L,K},"ClosedLoopSystem"];
initialConditions={0,1,0,0};                (* initialize system at {0,1} for δ-kick and estimator at {0,0}, no other direct inputs *)
{θ1,θd1,θh1,θhd1}=StateResponse[{closedLoop,initialConditions},0,t];   
p1a=Plot[{θ1,θh1},{t,0,tmax},PlotRange->All,ImageSize->250];
p2a=Plot[{θd1,θhd1},{t,0,tmax},PlotRange->All,ImageSize->250];
Grid[{{p1a,p2a}},Spacings->{4}]

{closedLoop,G1ss}  (* dynamics are the same; what is considered an input / output is slightly different  *)

Miscellaneous quantities for the problem:

TraditionalForm/@{θ,θd,θh,θhd}

MatrixForm /@{Ac,K,L,A}

Kob = K.Inverse[s(({{1, 0}, {0, 1}}))-A2].L // Simplify;
KobTF=TransferFunctionModel[Kob,s]