Using Mathematica to Design a Linear Observer Circuit for a Hyperchaotic Oscillator

Damon A. Miller			
Western Michigan University	    
d.a.miller@ieee.org

Giuseppe Grassi
University of Lecce

Introduction

Observers estimate target system states and thus are widely used in control and communication systems. This abstract presents the application of Mathematica to the design of a linear observer circuit for a hyperchaotic oscillator. Mathematica is used to simulate the hyperchaotic oscillator, determine state feedback gains for pole (eigenvalue) placement, investigate the sensitivity of system stability to these gains, verify observer performance, investigate noise effects, and insure that control signals are within the dynamic range of circuit elements. This abstract provides all needed code and other details of using Mathematica not presented in the companion paper [1]. That paper describes more general theoretical results and detailed experimental data for the observer designed here.

Hyperchaotic Oscillator

The following packages and definitions are required before proceeding:

<< LinearAlgebra`MatrixManipulation`; << 
      Graphics`Graphics`;
mH = 0.001; uF = 10^(-6.0); Kohms = 1000; mV = 0.001; 
    mA = 0.001; ms = 0.001;

The hyperchaotic oscillator of interest here was proposed in [2] and is shown in Figure 1 . This circuit is described by a set of differential equations

(({{(Subscript[(v'), C1])}, {(Subscript[(v'), C2])}, {i' L1}, {i' L2}}))=(({{(1/(R1 C1)), 0, ((-1)/C1), ((-1)/C1)}, {0, 0, 0, (1/C2)}, {(1/L1), 0, ((-RL1)/L1), 0}, {(1/L2), ((-1)/L2), 0, ((-RL2)/L2)}}))(({{(Subscript[v, C1])}, {(Subscript[v, C2])}, {(Subscript[i, L1])}, {(Subscript[i, L2])}}))+(({{0}, {((-1)/(R2 C2))}, {0}, {0}}))f[vC2]

                = A (({{(Subscript[v, C1])}, {(Subscript[v, C2])}, {(Subscript[i, L1])}, {(Subscript[i, L2])}}))+B f[vC2] = A x + B f[vC2]

where f[(Subscript[v, C2])]=R2 iD. Note that A is a 4x4 matrix of component values and x is the state vector to be estimated.

Show[Import["4DOsc.tif"]];

Figure 1. Four dimensional hyperchaotic oscillator [2]

Diode Model

The most straightforward approach to simulating this circuit is to numerically solve the associated differential equations (1).  The relationship between the diode current iD and the node voltage vC2 is

vC2 = R2 iD + vD[iD]

where vD is the forward diode voltage drop. Explicitly solving for iD[vC2] is not possible as vD[iD] is a transcendental function. In order to approximate this function an interpolation over numerically generated solution points is used. First, the component values for the drive system (the hyperchaotic oscillator), the observer system (to be discussed), and the "ideal" values of the components are set.

actualdrivevals={L1->98.3 mH,L2->32.6mH,C1->0.0980uF,C2->0.0327uF,
  R1->1.5594 Kohms,R2->50.5,RL1->68.5,RL2->32.3};
actualobsvals={L1->98.8 mH,L2->33.1mH,C1->0.0976uF,C2->0.0326uF,
  R1->1.5740 Kohms,R2->50.5,RL1->70.3,RL2->32.4};
standardvals={L1->100.0 mH,L2->33.0mH,C1->0.1uF,C2->0.033uF,
  R1->1.6 Kohms,R2->51,RL1->69.15,RL2->32.3};

This section of code generates an approximate f[(Subscript[v, C2])].

iD[vD_,Is_,n_, VT_]:=Is(Exp[vD/(n VT)]-1);
Off[Solve::ifun];Off[InverseFunction::ifun];
solA=Solve[iD[0.7,Is,2,25.2 mV]==4.72 mA, Is][[1]];
vD[iD_, Is_, n_, VT_]:=n VT Log[iD/Is+1];
iD[vC2_]:=iD /. Solve[(R2 /. actualdrivevals) iD  + 
           vD[iD, Is /. solA, 2, 25.2 mV]==vC2,iD][[1]];
tmp=Table[{vC2,iD[vC2] (R2 /. actualdrivevals)},{vC2,0,1,0.01}];
zeropad=Table[{vC2,0},{vC2,-2.5,-0.1,0.1}];
iDTable=Join[zeropad,tmp];
f=Interpolation[iDTable,InterpolationOrder->1];

Simulation

Using the approximate f[(Subscript[v, C2])] system (1) can be numerically solved; these results will be needed to simulate the observer. Note that time is normalized to 0.16 ms. Related plots will be provided later in the abstract.

Off[General::spell1];
A={{1/(R1 C1), 0, -1/C1, -1/C1},
   {0,0,0,1/C2},
   {1/L1,0,-RL1/L1,0},
   {1/L2,-1/L2,0,-RL2/L2}};
B={0,-1/(R2 C2),0,0};
tnorm=(R1 C1 /. standardvals);
Anormsym=tnorm A;Bnormsym=tnorm B;
Anorm=Anormsym /. actualdrivevals;
Bnorm=Bnormsym /. actualdrivevals;
x={vC1[t],vC2[t],iL1[t],iL2[t]};
tmp=Anorm . x + Bnorm f[vC2[t]];vars=x;
eqns={vC1'[t]==Part[tmp,1],vC2'[t]==Part[tmp,2],
      iL1'[t]==Part[tmp,3], iL2'[t]==Part[tmp,4],
      vC1[0]==-0.05,vC2[0]==-0.6,iL1[0]==-0.0002,
      iL2[0]==0.0002};
actualdrivesol=NDSolve[eqns,vars,{t,0,170},MaxSteps->25000];

Linear Observer

The linear observer circuit is shown in Figure 2.

Show[Import["LinObs.tif"]];

Figure 2. Linear observer for the hyperchaotic oscillator of Figure 1.

Design via Pole Placement

If the hyperchaotic oscillator supplies the scalar signal y = f[vC2] + K x, where K is a vector of state feedback gains, and the observer supplies the signal yhat = K xhat, where xhat contains the analogous observer states,  the error between the two systems is linear and described by [1]

e' = (x-xhat )' = (A - B K) e .

The observer error e can be asymptotically forced to zero by an appropriate choice of K since K sets the error system pole (eigenvalue) locations. These pole locations are chosen to provide an ITAE response [3] which corresponds to the characteristic equation:

Chop[Expand[w0^4 (s/w0 + 0.4240 + 1.2630 I)(s/w0 + 0.4240 - 1.2630 I)(s/w0+0.6260 + 0.4141 I)(s/w0 + 0.6260 - 0.4141 I)]]

The next section of code generates a table of gains for the normalized system for various values of w0.

Anorm= Anormsym /. standardvals;Bnorm=Bnormsym /. standardvals;
ITAE[s_,w0_]:=
  s^4 + 2.1 s^3 Sqrt[w0] + 3.4s^2 Sqrt[w0]^2 + 
    2.7 s Sqrt[w0]^3+Sqrt[w0]^4;
gains[w0_] :=
   (K=Array[k,4];
    desiredpoly=ITAE[s,w0];                        
    charpoly = Chop[Det[s IdentityMatrix[4]-
      (Anorm-Outer[Times,Bnorm,K])]];
    eqns={Coefficient[desiredpoly,s]==Coefficient[charpoly,s],                                     
          Coefficient[desiredpoly,s^2]==Coefficient[charpoly,s^2],                                     
          Coefficient[desiredpoly,s^3]==Coefficient[charpoly,s^3],                                     
          Part[desiredpoly,1]==
          Part[charpoly,1]+Part[charpoly,6]+Part[charpoly,8]+
          Part[charpoly,12]+Part[charpoly,13]};
          Ksol=Solve[eqns,K][[1]]);
 Table[gains[w0],{w0,10,50,10}] // MatrixForm

The gains for w0=30 rad/sec are chosen and the resulting eigenvalues are checked:

idealKsol=gains[30]
Eigenvalues[Anorm-Outer[Times,Bnorm,K /.idealKsol]]

Realistic gains are used for subsequent calculations yielding these normalized and unnormalized eigenvalues:

actualKsol={k[1]->-0.555,k[2]->-0.129,k[3]->-502,k[4]->168};
Eigenvalues[Anorm-Outer[Times,Bnorm,K /.actualKsol]]
Eigenvalues[Anorm-Outer[Times,Bnorm,K /.actualKsol]] / (0.16 ms)
Eigenvalues[(A /. standardvals)-Outer[Times,(B /. standardvals),K /.  actualKsol]]

Finally, error system pole locations for 500 gains with a uniform 5% error are checked. Poles remain in the left half plane and thus the error system remains stable for this level of gain error.

k1l=0.95 k[1] /. actualKsol; k1u=1.05 k[1] /. actualKsol;
k2l=0.95 k[2] /. actualKsol; k2u=1.05k[2] /. actualKsol;
k3l=0.95k[3] /.actualKsol; k3u=1.05 k[3] /. actualKsol;
k4l=0.95 k[4] /.actualKsol;k4u=1.05 k[4] /. actualKsol;
For[{i=1;xypts={};SeedRandom[5]},i<=500,i++,
	{randK={k[1]->Random[Real,{k1l,k1u}],
	        k[2]->Random[Real,{k2l,k2u}],
			k[3]->Random[Real,{k3l,k3u}],
			k[4]->Random[Real,{k4l,k4u}]};
	 eigenvals=
      Eigenvalues[(A /. standardvals) - 
          Outer[Times,(B /. standardvals),(K/. randK)]];
   xypts=Insert[xypts,{Re[Part[eigenvals,1]],Im[Part[eigenvals,1]]},1];
   xypts=Insert[xypts,{Re[Part[eigenvals,2]],Im[Part[eigenvals,2]]},1];
   xypts=Insert[xypts,{Re[Part[eigenvals,3]],Im[Part[eigenvals,3]]},1];       
   xypts=Insert[xypts,{Re[Part[eigenvals,4]],Im[Part[eigenvals,4]]},1];}
 ]; 
 ListPlot[xypts, Frame->True, AxesLabel->{"Re[s]","Im[s]"}, 
         AxesOrigin->{0,0}, 
         PlotRange->{{-30000,1000},{-50000,50000}},
         FrameTicks->
         {{-10000,-20000,1000},{-40000,-20000,0,20000,40000},{},{}},
         FormatType->OutputForm, Axes->True];

Figure 3. Eigenvalues for 500 random feedback gain vectors having a uniform 5% error

This completes the observer design.

Simulation

This section verifies the observer performance for two cases: a noise-free  and a noisy y[t]. Both simulations use actual component values.

Noise Free y[t]

The resulting synchronizing signal is not too large for circuit implementation:

y1=Head[(vC1[t] /. actualdrivesol)[[1]]];
y2=Head[(vC2[t] /. actualdrivesol)[[1]]];
y3=Head[(iL1[t] /. actualdrivesol)[[1]]];
y4=Head[(iL2[t] /. actualdrivesol)[[1]]];
y[t_]:=((k[1]/.actualKsol) y1[t] +
		(k[2]/.actualKsol) y2[t]+ 
		(k[3]/.actualKsol) y3[t]+
		(k[4] /. actualKsol)y4[t] +	
	     f[y2[t]]);
Plot[y[t],{t,50,80}, AxesLabel->{"t","y"}];

Figure 4. The synchronizing signal

This code simulates the observer and plots the oscillator (blue) and observer (red) vC2 state. Note that the observer closely follows the oscillator.

Anorm=Anormsym /. actualobsvals;
Bnorm=Bnormsym /. actualobsvals;
xhat={vC1hat[t],vC2hat[t],iL1hat[t],iL2hat[t]};
yhat[t_]:=(k[1] /. actualKsol) vC1hat[t] + (k[2] /. actualKsol) vC2hat[t] + (k[3] /. actualKsol) iL1hat[t] +(k[4]/.actualKsol)iL2hat[t];
tmp = Anorm . xhat + Bnorm (y[t]-yhat[t]);
vars=xhat;
eqns={vC1hat'[t]==Part[tmp,1],vC2hat'[t]==Part[tmp,2],
iL1hat'[t]==Part[tmp,3],iL2hat'[t]==Part[tmp,4],
vC1hat[0]==-0.01,vC2hat[0]==0,iL1hat[0]==0,iL2hat[0]==0};
actualobssol=NDSolve[eqns,vars,{t,0,170},MaxSteps->120000];
Plot[
  {Evaluate[vC2hat[t]/.actualobssol],Evaluate[vC2[t]/.actualdrivesol]},
  {t,50,80},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}},
  PlotRange->All,AxesLabel->{"t","vC2"}];

Figure 5. Hyperchaotic oscillator and observer waveforms for state vC2

Noisy y[t]

The effects of a noisy y[t] are investigated by adding this noise signal to y[t]:

mean = 0; variance = (20  mV)^2; 
    SeedRandom[10];
noisedata = 
        Table[{t, 
            mean + Sqrt[3.0/0.5^2] (Random[] - 0.5) Sqrt[variance]}, {t, 
            0, 200, 0.05}];
noise = Interpolation[noisedata];
Plot[noise[t], {t, 50, 80}, AxesLabel -> {"t", "noise"}, 
        PlotRange -> All];

Figure 6. Noise signal

This clearly degrades the observer performance:

Anorm=Anormsym /. actualobsvals;
Bnorm=Bnormsym /. actualobsvals;
xhat={vC1hat[t],vC2hat[t],iL1hat[t],iL2hat[t]};
yhat[t_]:=(k[1] /. actualKsol) vC1hat[t] + (k[2] /. actualKsol) vC2hat[t] + (k[3] /. actualKsol) iL1hat[t] +(k[4]/.actualKsol)iL2hat[t];
tmp = Anorm . xhat + Bnorm (y[t]-yhat[t]+noise[t]);
vars=xhat;
eqns={vC1hat'[t]==Part[tmp,1],vC2hat'[t]==Part[tmp,2],
iL1hat'[t]==Part[tmp,3],iL2hat'[t]==Part[tmp,4],
vC1hat[0]==-0.01,vC2hat[0]==0,iL1hat[0]==0,iL2hat[0]==0};
actualobssol=NDSolve[eqns,vars,{t,0,170},MaxSteps->120000];
Plot[
  {Evaluate[vC2hat[t]/.actualobssol],Evaluate[vC2[t]/.actualdrivesol]},
  {t,50,80},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}},
  PlotRange->All,AxesLabel->{"t","vC2"}];

Figure 7. Hyperchaotic oscillator and observer waveforms for state vC2

Conclusion

This abstract has described how Mathematica was used to design a linear observer circuit for a hyperchaotic system. Consult [1] for further details including a quantitative evaluation of simulated and experimental observer performance.

References

[1]  D. A. Miller and G. Grassi, "Experimental realization of observer-based  hyperchaos synchronization," IEEE Transactions on Circuits and Systems-I,  vol. 48, pp. 366-374,  March 2001.

[2] A. Tamasevicius, A. Namajunas, and A. Cenys, "Simple 4D chaotic oscillator," Electronics Letters, vol. 32, no. 11, pp. 957-958, 1996. 

[3] G. F. Franklin, J. D. Powell, and A. Emami-Naeini, Feedback Control of Dynamic Systems. Reading, MA: Addison-Wesley, 1986.