SINGULAR PERTURBATION
     NCAlgebra NOTEBOOK
           follows CDC99
         by  Dell   and Bill

The system we study and its control (Sec 2)

(({{(dx/dt)}, {ϵ  (dz/dt)}})) = (({{(Subscript[A, 11]), (Subscript[A, 12])}, {(Subscript[A, 21]), (Subscript[A, 22])}}))(({{x}, {z}})) + (({{(Subscript[B, 1])}, {(Subscript[B, 2])}})) u
y = (({{(Subscript[M, 1]), (Subscript[M, 2])}}))(({{x}, {z}}))
	Minimize   J  =   
          (∫_0^∞)   (y^T)  y  +  
            (u^T) R u     dt
\(where\ \ R\ is\ a\ positive\ definite\ weight\ matrix . 
It\ is\ well\ 
        known\ that\ the\ solution\ to\ this\ problem\)is of the form 
\(u\^*\) =  - (R^(-1))  (B^T)  K(ϵ) (({{x}, {z}}))
where  K is the solution to the ARE
        K A   +  (A^T) K  -  K B  (R^(-1))   (B^T) K   + 
             (M^T)  M  =  0
with 
A =    (({{(Subscript[A, 11]), (Subscript[A, 12])}, {((Subscript[A, 21])/ϵ), ((Subscript[A, 22])/ϵ)}})),  B =   (({{(Subscript[B, 1])}, {((Subscript[B, 2])/ϵ)}})),  
           and   M = (({{(Subscript[M, 1]), (Subscript[M, 2])}}))
K  =  
          (Subscript[K, 0])  +    ϵ (Subscript[K, 1])   +  (ϵ^2)  (Subscript[K, 2])  + 
                  ...  .

Format K compatibly with the partitioning of A using the notation
of [KKO86].
		K = (∑_(i = 0)^∞) (ϵ^i) (({{(Subscript[k, 1  i]), ϵ (Subscript[k, 2  i])}, {ϵ (Subscript[k, 2  i]^T), ϵ (Subscript[k, 3  i])}}))

Computer Algebra vs the problem (Sec3)

The ZEROth order term of the Riccati equation (Sec3.1)

Computer Algebra Finds the Basic Equation (Sec 3.1.1)

We begin by  setting the matrices in our problem to 
be noncommutative.

SetNonCommutative[
	k10,k20,k30,A11, 
	A12,A21,A22,B1,B2,M1,M2,R
	];

Next we create the first term of the power series 
      
         K  =   (Subscript[K, 0])  + ϵ (Subscript[K, 1]) +  (ϵ^2)  (Subscript[K, 2])  +   ...  . 
        

It contains an ϵ.

K0 = {{k10, ep**k20},
     {ep**tp[k20], ep**k30 }};
     
K0//MatrixForm

Next we define the block matrices present in our problem.

A  = {{A11, A12},
      {(1/ep)**A21,(1/ep)**A22}};
B  = {{B1},{(1/ep)**B2}};
Mm = {{M1,M2}}; 

A//MatrixForm 
B//MatrixForm
Mm//MatrixForm

The following Mathematica input creates a function which 
makes Riccatis and sets G to be the optimal control.

Riccati[K_] := MatMult[K,A] + 
	MatMult[tpMat[A],K ] -
	MatMult[K,B,{{Inv[R]}},tpMat[B],K] +
	MatMult[tpMat[Mm],Mm]

G = - MatMult[{{Inv[R]}},tpMat[B],K0]

Extract the zero-th order terms in ϵ 
from the upper left entry of the perturbed
Riccati  to create a polynomial.

Ep10 = NCTermsOfDegree[ 
			Riccati[K0][[1,1]],{ep},{0} ]

Extract the zero-th order terms in ϵ from the upper right entry of the perturbed Riccati  to create a polynomial.

Ep20 = NCTermsOfDegree[ 
			Riccati[K0][[1,2]],{ep},{0} ]

NCCollect[ NCExpand[Ep20], {k20}]

Notice that by assuming the invertibility of  
                               (A22-B2**Inv[R]**tp[B2]**k30)
 in the above equation we can solve for k20.

         k20 = ( k10**A12  +tp[A21]**k30+
               tp[M1]**M2-k10**B1**Inv[R]**tp[B2]**k30 ) **   
               Inv[ A22-B2**Inv[R]**tp[B2]**k30 ]

This observation will be used below.

Next, extract the zero-th order terms in ϵ  from the 
lower right entry of the perturbed Riccati  to create a polynomial.

Ep30 = NCTermsOfDegree[ 
			Riccati[K0][[2,2]],{ep},{0} ]

Notice that the above equation is a Riccati in k30 and can be used to
determine this matrix.

The matrix which is still unknown is k10.  Solving for it by hand is a bit tricky, see Section 3.4 [KKO86].
Solving for k10 by computer is what  we turn to after introducing  other basic though rinky dink algebraic identities.

Other Algebraic Identities Known to Hold (Sec 3.1.3)

Abbreviations which are standard conventions in the field [KKO86]:

SetNonCommutative[G1,G2,N0,M0,A0,B0,R0];
Abbreviations = {
	G1 == G[[1,1]]/.ep->0,
	G2 == G[[1,2]]/.ep->0,
	N0 == - M2**Inv[A22]**B2,
	M0 == M1 - M2**Inv[A22]**A21,
	A0 == A11 - A12 ** Inv[A22]**A21,
	B0 == B1 - A12**Inv[A22]**B2,
	R0 == R + tp[N0]**N0
	};

Also some of our matrices are known to be self adjoint.

SelfAdjoint = {
	k10 == tp[k10],
	k30 == tp[k30],
	R == tp[R],
	R0 == tp[R0],
	Inv[R] == tp[Inv[R]],
	Inv[R0] == tp[Inv[R0]]
	};

It is standard that A0, A11, and A22 are invertible.
Also, the fast  sub-system is (assumed ) strictly stable which makes
                   A22 + B2**G2  invertible.  
Recall  the 
              invertibility of    A22 - B2**Inv[R]**tp[B2]**k30 
 follows from our solution for k20.
We make the resulting relations with the next command.

inverses = NCMakeRelations[{Inv,
	R,R0,A0,A11,A22,(A22 + B2**G2), 
	(A22 - B2**Inv[R]**tp[B2]**k30) }];

Solving for k10 using a  Grobner Basis.

Next we combine all the relations created above and add in 
their tranposes.

Relations = Union[{Ep10,Ep20,Ep30},
 	Abbreviations,
	SelfAdjoint, inverses ];

AllRelations = NCAddTp[ Relations ]

Now that we have created all equations which need to be solved 
simultaneously to find k10,  we invoke the NC Groebner Basis
program NCProccess to solve for  k10. First we set the order reflecting our priorities on which variables should be eliminated.

NCAutomaticOrder[ {{N0,M0,R0,A0,B0}, 
		{k10},{B1,B2,M1,M2,R,
		A11,A12,A21,A22,
		Inv[A22-B2**Inv[R]**tp[B2]**k30],
	    tp[Inv[A22-B2**Inv[R]**tp[B2]**k30]]},
	    {k30},{k20},{G1,G2}},  AllRelations];

Finally we make the call to NCProcess[ ]  which will 
create a sorted Groebner basis from the relations in AllRelations
and hopefully give us an answer for k10.  The Mathematica answer
can be found below and the TeX answer can be found on the 
main Singular Perturbation page.

answer = NCProcess[AllRelations,4,
	"FindK10",SBByCat->False  ,
	DegreeCap->14, DegreeSumCap->17] ;

answer