Miscellaneous explorations of model reduction.  
Section 6.5.2 (Balanced coordinates)  Figure 6.7

• revised, Oct. 27, 2023:
  – use SystemsModelDelete to simplify  model-reduction commands
  – small changes to make code or output prettier

Clear["Global`*"];
SetOptions[Plot,PlotRange -> All];

A1=(({{-1, 0}, {0, -2}})); B1=(({{1}, {1}})); C1=(({{1, 2}}));
Gss1=StateSpaceModel[{A1,B1,C1}]

TransferFunctionModel[Gss1,s]

Apart[({{((4+3 s)/(2+3 s+(s^2)))}})]

Gsstr1=SystemsModelDelete[Gss1,None,None,2]        (* delete the second mode *)

{atr,btr,ctr,dtr}=Normal[Gsstr1];            (* get the explicit matrices *)

dres=DiagonalMatrix[-C1[[1,2]]/A1[[2,2]]*B1[[2]]];
Gssres1=StateSpaceModel[{atr,btr,ctr,dres}]

p1=ControllabilityGramian[Gss1];q1=ObservabilityGramian[Gss1];TraditionalForm/@{p1,q1}

{tbal1,Gssbal1}=InternallyBalancedDecomposition[Gss1]//FullSimplify;
{abal1,bbal1,cbal1,dbal1}=Normal[Gssbal1];
TraditionalForm /@ {tbal1,Gssbal1}

Gssbaltr1=SystemsModelDelete[Gssbal1,None,None,2]
{atr1,btr1,ctr1,dtr1}=Normal[Gssbaltr1];

ares1=DiagonalMatrix[{abal1[[1,1]]-abal1[[1,2]]/abal1[[2,2]]*abal1[[2,1]]}]//FullSimplify;
bres1=DiagonalMatrix[{bbal1[[1,1]]-abal1[[1,2]]/abal1[[2,2]]*bbal1[[2,1]]}]//FullSimplify;
cres1=DiagonalMatrix[{cbal1[[1,1]]-cbal1[[1,2]]/abal1[[2,2]]*abal1[[2,1]]}]//FullSimplify;
dres1=DiagonalMatrix[{-cbal1[[1,2]]/abal1[[2,2]]*bbal1[[2,1]]}]//FullSimplify;
Gssbalres1=StateSpaceModel[{ares1,bres1,cres1,dres1}]

TransferFunctionModel /@ {Gss1,Gsstr1,Gssres1,Gssbaltr1,Gssbalres1}//Simplify

{BodePlot[{Gss1,Gsstr1,Gssres1}],BodePlot[{Gss1,Gssbaltr1,Gssbalres1}]}

y1=OutputResponse[Gss1,UnitStep[t],t];
y2=OutputResponse[Gsstr1,UnitStep[t],t];
y3=OutputResponse[Gssres1,UnitStep[t],t];
y4=OutputResponse[Gssbaltr1,UnitStep[t],t];
y5=OutputResponse[Gssbalres1,UnitStep[t],t];
{Plot[{y1,y2,y3},{t,0,5}],Plot[{y1,y4,y5},{t,0,5}]}

p1bal=ControllabilityGramian[Gssbal1];
q1bal=ObservabilityGramian[Gssbal1];
TraditionalForm /@ {p1bal,q1bal,p1bal//N}

Eigenvalues[Gssbal1[[1,1]]]

FullSimplify /@{y1,y2,y3,y4,y5}//TraditionalForm

Look at 1d optimal Hankel approximant (computed in Scilab)

ah={{-(Sqrt[2])}}; bh={{1+(1/(Sqrt[2]))}}; ch={{(2/3)(1+(Sqrt[2]))}};dh={{(1/2)-((Sqrt[2])/3)}};
Gssh=StateSpaceModel[{ah,bh,ch,dh}]

y6=OutputResponse[Gssh,UnitStep[t],t];
LogLogPlot[{y1,y6},{t,0.01,100}]

y6//FullSimplify