APPENDIX 5

Crane

Imported Packages

<< Miscellaneous`PhysicalConstants`

Needs["Statistics`ContinuousDistributions`"]

Naming Conventions

In order to simplify the code following naming conventions were adopted: 
	- coordinates of vectors have the form x[i] or (Subscript[x, i]) on the input. The later form is used on output. Whole vectors are denoted as x[]
	- the letters z[i], y[i] can be exceptionally used in the same way
	- time is always denoted as t
	- integration constant and a integration functions for PDEs are denoted by C1, C2 ...
	- vectorfields are usually called f,g
	- distributions D
	- coordinate transformations S,T,W

Before any vectors are used, they must be defined. At the same time we establish equivalence between x[i] or (Subscript[x, i]).

DefineVector[x_, n_]  :=  
      Module[{i}, x[] = Array[x, n]; Do[x[i] = (Subscript[x, i]),  {i, n}];  
        MaxOrder  =  n]

The DSolve[] function does not work correctly with indexed variables like x[1], x[2]. Therefore we need a pair of textual transformations to replace any x[1] with a$$1 before calling the DSolve and to replace bakwards a$$1 with x[1] after calling DSolve. We will transform integration constants C[i] to Cxx[i] to make them unique in case more than one DSolve result will be combined into one expression. Unfortunatelly, we encountered problems with the option DSolveConstant, thats why the operation must be made by the hand.

We assume that only the x[] vector should be transformed. Maybe we will add options in the future.

\(\(DSolvePreprocess[chain_] := 
      Module[{i}, \ 	
		
        chain\  /. \ 
          Table[x[i]\  -> \ Symbol["a$$"<>ToString[i]], {i, MaxOrder}]]	
    \)\)

\(DSolvePostprocess[chain_] := 
      Module[{i}, 
		
        \((chain\  /. \ 
              Table[Symbol["a$$"<>ToString[i]]\  -> \ x[i], {i, 
                  MaxOrder}])\)
		 /. 
          C[a_]\  :> \ Symbol["C"<>ToString[a]]
	]\)

\(MyDSolve[eqn_, y_, x_]\  := \ 
      Module[{i}, 
		
        DSolvePostprocess[
          DySolve[
				DSolvePreprocess[eqn], DSolvePreprocess[y], \ 
            DSolvePreprocess[x]
			]]]\)

Initial conditions

InitialConditions[n_]  :=  Table[(Subscript[x, i])[0] ==  (Subscript[x0, i]), {i, n}]

Basic Differential Operators

Gradient and  Jacobian.

DScalarByVector[sc_, vec_] :=  Map[D[sc, #]&, vec]

DVectorByVector[vec_, byvec_]  :=  
      Map[DScalarByVector[#, byvec]&, vec]

Lie derivative and Lie bracket

LieBracket[f_, g_, X_]  :=  
      Simplify[Dot[DVectorByVector[g, X], f] - Dot[DVectorByVector[f, X], g]]

LieDerivative[f_, h_,  X_]  :=  
      Simplify[Dot[f,  DScalarByVector[h, X]]]

Multiple Lie derivative and Lie bracket

LieBracketNN[f_, g_, X_, n_]  :=  
      If[n  <=  0, g, LieBracket[f, LieBracketNN[f, g, X, n - 1], X]]

LieDerivativeNN[f_, h_, X_, n_]  :=  
      If[n  <=  0, h, LieDerivative[f, LieDerivativeNN[f, h, X, n - 1], X]]

Flow

\(Flow[f_]\  := \ 
      Module[{i, n, \ f1, \ eq}, \ 
		n\  = \ Length[f]; 
		
        f1\  = \ f\  /. \ x\_a_\  -> \ x\_a[t]; 
		
        eq\  = \ Join[Table[f1[\([i]\)]\  == \(x\_i'\)[t], {i, n}], 
            InitialConditions[n]]; 
		
        MyDSolve[eq, Table[x\_i[t], {i, n}], t]]\)

Vectorfield transformation

TStar[f_, S_]  :=  (DVectorByVector[S, x[]]  .  f)

Crane

State equations

\(\(
\(eq1\  = \ 
      \((ml\  + \ mc)\)\ dx4\  + \ 
          ml\ l\ \((dx2\ Cos[x[1]]\  - \ x[2] Sin[x[1]])\)\  == \ u; \)\)\)

eq2  =  
      (ml dx4 Cos[x[1]]  +  ml l dx2  ==  
          -ml Gravity Sin[x[1]]);

eq3  =  dx2  /.  Solve [eq2  /.  Solve[eq1,  dx4],  dx2];

dxx2  =  Simplify[eq3[[1]]];

g2  =  Coefficient[dxx2, u]

f2  =  Simplify[dxx2  -  u g2]

Complete state space model

DefineVector[x, 2];

f = {x[2], f2};  g = { 0,  g2}; 
    σ = {0,  F Cos[x[1]] / (ml l)};

Correcting term

dσ  =  DVectorByVector[σ, x[]]

corr  =  dσ  .  σ

Check if the distribution is nonsingular (involutivness is assured for n <= 2).

ad = LieBracket[f, g, x[]]

NullSpace[{ad, g}]

Now select the output function and check the relative degree

h  =  x[1];

LieDerivative[g, h, x[]]

LieDerivative[g, LieDerivative[f, h, x[]], x[]]

Compute transformation

λ = Integrate[F/(ml l σ[[2]]), x[1]]

LieDerivative[g, λ, x[]]

T2  =  LieDerivative[f, λ, x[]]

b  =  LieDerivative[g, T2, x[]]

a  =  LieDerivative[f,  T2,  x[]]

T = {λ, T2}

Vectorfields fo the system in the new coordinates z

\(\(f1\  = \ Simplify[TStar[f, T]]
\)\)

g1  =  Simplify[TStar[g, T]]

σ1  =  Simplify[TStar[σ, T]]

Compute the feedback

β = 1/b

α = -a/b

System with the feedback

\(f\&^ = Simplify[f1 + \ α\ g1]\)

\(g\&^ = g1\ β\)

Numerical experiments

values  = 
       { ml  ->  1000,  mc  ->  10000,  l  ->  30,  F -> 250,  
        Gravity  ->  9.81};

{TT,  ff, gg, σσ,  αα,  ββ}
          =  {T,  f, g,  σ,  α,  β}  /.  values;

Plot[TT[[1]]  /.  x[1]  ->  t, {t, -Pi/2,  Pi/2}]

Plot[TT[[2]]/x[2]  /.  x[1]  ->  t, {t, -Pi/2,  Pi/2}]

Plot[ββ /.  x[1]  ->  t, {t, -Pi/2, Pi/2}]

Period of equivalent mathematical pendulum

Tmath  =  Pi Sqrt[l / Gravity]

N[Tmath  /.  values]

Simulation

\(CraneExperiment[vt_, wt_, \ dt_, \ x0_, \ k_]\  := \ 
      Module[{n, i, \ va}, 
		n\  = \ Length[vt]; 
		
        xt\  = \ Table[0, {i, n}]; 			
		zt\  = \ Table[0, {i, n}]; 

        \ \ \ \ ut\  = \ Table[0, {i, \ n}]; 
		xt[\([1]\)]\  = \ x0; 

        		For[i = 1, i < n, \(i++\), 
			\ 
          va\  = \ {x[1]\  -> \ \(xt[\([i]\)]\)[\([1]\)], \ 
              x[2]\  -> \ \(xt[\([i]\)]\)[\([2]\)]}; \ 
		\ 	
          zt[\([i]\)]\  = \ N[TT\  /. \ va]; 
\ \ \ \ \ \ \ 
          ut[\([i]\)]\  = \ \ 
            N[αα\ \  + 
                  ββ\ \ 
                    \((vt[\([i]\)]\  + \ 
                        k[\([1]\)] \(zt[\([i]\)]\)[\([1]\)]\  + \ 
                        k[\([2]\)]\ \ \(zt[\([i]\)]\)[\([2]\)])\)\ \  /. \ 
                va]; 				\ 	\ \ 
          dx\  = \ N[
              \((dt\ \((\ 
                      ff\  + \ gg\ ut[\([i]\)]\  + \ 
                        σσ\ wt[\([i]\)])\) /. \ va)\)]; 
			
          xt[\([i + 1]\)]\  = \ xt[\([i]\)] + dx; 
		]; 
		
        zt[\([n]\)]\  = \ zt[\([n - 1]\)]; 
		
        ut[\([n]\)]\  = \ ut[\([n - 1]\)]; 
	]\)

\(\(DoExperiment[n_, vt_, wt_, k_]\  := \ 
      Module[\ {i}, 
	dt\  = \ 0.01; 
	
        tt\  = \ Table[dt\ \((i - 1)\), \ {i, n}]; 	
\ \ 
        CraneExperiment[vt, wt, 1/100, \ {0, 0}, \ k]; 
\ \ 
        x1l\  = \ \(Transpose[xt]\)[\([1]\)]; 
	
        x2l\  = \ \(Transpose[xt]\)[\([2]\)]; 
\ \ 
        z1l\  = \ \(Transpose[zt]\)[\([1]\)]; 
	
        z2l\  = \ \(Transpose[zt]\)[\([2]\)]; 
\ \ 
        pz1 = \ ListPlot[
            Transpose[{tt, z1l}], {PlotLabel\  -> \ "z1", 
              DisplayFunction -> Identity}]; 
\ \ 
        pz2\  = \ 
          ListPlot[
            Transpose[{tt, z2l}], {PlotLabel\  -> \ "z2", 
              DisplayFunction -> Identity}]; 
\ \ 
        px1 = \ ListPlot[
            Transpose[{tt, x1l}], {PlotLabel\  -> \ "x1", 
              DisplayFunction -> Identity}]; 
\ \ 
        px2\  = \ 
          ListPlot[
            Transpose[{tt, x2l}], {PlotLabel\  -> \ "x2", 
              DisplayFunction -> Identity}]; 
	
        pu\  = \ ListPlot[
            Transpose[{tt, ut}], {PlotLabel\  -> \ "u", 
              DisplayFunction -> Identity}]; 
	
        pw\  = \ ListPlot[
            Transpose[{tt, wt}], {PlotLabel\  -> \ "w", 
              DisplayFunction -> Identity}]; 
        Show\ [GraphicsArray[{{pz1, \ pz2}, {px1, \ px2}, \ {pu, \ pw}}]]; 

        	]; \)\)

Deterministic experiment, no noise, ramp input, no control

DoExperiment[150, Table[1, {i, 150}],  Table[0,  {i, 150}],  {0, 0}]

Stochastic experiment,  noise, pulse input, no control

pulse  =  Table[0, {i, 150}];  pulse[[1]] = 100;

DoExperiment[150, pulse,  Table[0,  {i, 150}],  {-10, -5}]

Stochastic experiment, noise, no feedback

DoExperiment[150, Table[0, {i, 150}],  
      Table[Random[NormalDistribution[0, 1000]],  {i, 150}],  {0, 0}]

Stochastic experiment, noise,  feedback

DoExperiment[150, Table[0, {i, 150}],  
      Table[Random[NormalDistribution[0, 1000]],  {i, 150}],  {0, -1}]