CLOSE APPROACH TO CONSTRAINT SET:OPTIMAL CONTROL PROBLEM
				
by

Paul J. Scott

Introduction

This notebook performs the numerical calculations proposed in a paper with same title.  References below refer to this paper.  This is a shooting method in which one varies the inital values of the costate until the constraint set is reached.

There are six parts.  Part zero contains a Runge Kutta DE solver for a vector system of equations.  These were copied from a MathSource package by Roman E. Maeder. Part One sets of the optimal control problem in the same notation as the paper.   Part Two sets up variables for use in the Runge-Kutta solver.  Part Three is a function, (Module), for evaluating the right hand side of all DE's for use in the Runge-Kutta solver.  This function includes computing the control.  Part Four sets up initial values for the state y and z. Part Five is a shooting step where initial values of the costates are adjusted until the constraint set is reached.

Part Zero

This function performs one Runge Kutta step for vector right hand side function f with vector y for initial vector y0 and step size dt.

\(rkStep[f_, \ y_, \ y0_, \ dt_]\  := 
	
      Module[{\ k1, \ k2, \ k3, \ k4\ }, 
		
        k1\  = \ dt\ N[\ f\  /. \ Thread[y\  -> \ y0]\ ]; 
		
        k2\  = \ dt\ N[\ f\  /. \ Thread[y\  -> \ y0\  + \ k1/2]\ ]; 
		
        k3\  = \ dt\ N[\ f\  /. \ Thread[y\  -> \ y0\  + \ k2/2]\ ]; 
		
        k4\  = \ dt\ N[\ f\  /. \ Thread[y\  -> \ y0\  + \ k3]\ ]; 
		
        y0\  + \ \((k1\  + \ 2\ k2\  + \ 2\ k3\  + \ k4)\)/6
	]\)

This function performs a series of fixed step size Runge Kutta steps to generate a trajectory.

\(\(rk[f_List, \ y_List, \ y0_List, \ {t1_, \ dt_}]\  := 
	
      NestList[\ rkStep[f, \ y, \ #, \ N[dt]]&, \ N[y0], \ Round[N[t1/dt]]\ ]
          \  /; 
		Length[f]\  == \ Length[y]\  == \ Length[y0]
\)\)

BOTH FUNCTIONS BASED ON A PACKAGE BY ROMAN E. MAEDER IN MathSource.

Part One

We set up equations using the same notation as in the paper.

Declaring Variables

Number of y state variables (called n).

n  =  2

Declaration of y state variables.

yy  =  Table[ (Subscript[y, i]) ,  { i ,  1 ,  n } ]

Declaration of y- costate variables (called λ).

λλ  =  
      Table[  (Subscript[λ, i])  ,  { i , 1 ,  Length[ yy ] }  ]

Number of z state variables (called m).

m  =  2

Declaration of z-state variables.

zz  =  Table[  (Subscript[z, i])  ,  { i , 1 , m }  ]

Declaration of z-costate variables (called ξ).

ξξ  =  
      Table[  (Subscript[ξ, i])  ,  { i , 1 , Length[ zz ] }  ]

Declaration of control variables (called u).  See paper (1.3).

uu  =  Table[  (Subscript[u, i])  ,  { i , 1 ,  Length[ zz ] }  ]

The following "nuu" will be needed for numerical calculations.

nuu  =  
      Table[  (Subscript[nu, i])  ,  { i , 1 ,  Length[ zz ] }  ]

Declaration of K-T multipliers.

See paper (4.1)

llm  =  Table[  (Subscript[lm, i])  ,  { i , 1 , Length[ zz ] }  ]

The following "nllm" will be needed for numerical calculations.

nllm  =  
      Table[  (Subscript[nlm, i])  ,  { i , 1 , Length[ zz ] }  ]

See paper (4.2).

llp   =   
      DiagonalMatrix[  Table[  (Subscript[lp, i]) ,  {i, 1, Length[ zz ] } ] ]

TableForm[ llp ]

The following "nllp" will be needed for numerical calculations.

nllp   =   
      DiagonalMatrix[  Table[  (Subscript[nlp, i]) ,  {i, 1, Length[ zz ] } ] ]

Id   =   IdentityMatrix[ Length[ zz ] ]

TableForm[ Id ]

Defining the differential equation and constraint

Definition of rhs of ODE   F.  See paper (1.1)

Each component of F is defined.

(Subscript[f, 1])   =   -(Subscript[y, 1])   +   (Subscript[z, 1])

(Subscript[f, 2])   =   -(Subscript[y, 2])   +   (Subscript[z, 2])

The components are listed in vector F.

F  =  Table[  (Subscript[f, i])  ,  { i , 1 , Length[ yy ] }  ]

Definition of constraint function G.  See paper (1.2)

Each component of g is defined.

(Subscript[g, 1])  =   
      (Subscript[y, 1])^2   +   (Subscript[y, 2])^2   -  ((Subscript[z, 1])^2   +   (Subscript[z, 2])^2 )

(Subscript[g, 2])  =    (Subscript[y, 1])   -   (Subscript[z, 1])

The components are listed in a vector G.

G  =  Table[  (Subscript[g, i])  ,  { i , 1 , Length[ zz ] }  ]

The Jacobian Matrix (∂_y) G

DyG   =   
      Transpose[  
        Table[ (∂_((Subscript[y, i]))) G ,  {i, 1, Length[ yy ]} ]  ]

TableForm[ DyG ]

The Jacobian Matrix (∂_z) G

DzG   =   
      Transpose[  
        Table[ (∂_((Subscript[z, i]))) G ,  {i, 1, Length[ zz ]} ]  ]

TableForm[ DzG ]

We will need the following definition later.

Define of G2.  See paper (1.5)

G2  =  
      Table[  G[[ i ]]^2  ,  { i , 1 , Length[ G ] }  ]

Definition of Hamiltonian.  See paper (4.3)

\(H\  = \ \ 
      G\ \ .\ \((\ Id\  - \ llp\ )\)\ .\ G\ \  + \ 
			\ \ \ 
        \((\ λλ\ \ \  + \ \ Transpose[\ DyG\ ]\ .\ llm\ )\).\ 
          F\ \  + \ \ 
\ \ \ \ \ \ \ \ \ \ 
        \((\ ξξ\ \ \ \  + \ Transpose[\ DzG\ ]\ .\ llm\ )\).\ uu\)

Setting up the right hand side of state costate equations

"stcst"  consists  of  the  rhs of  all  differential  equations.  Its  first component 1 is for integrating the time equation  (∂_t) t  =  1.  The order of all the components is t,y, z, λ, ξ.

First we include the rhs of (1.1) and (1.3).

stcst  =  Flatten[ { 1 ,  F  ,  uu } ]

Next we set up and append the rhs of system  (3.2)

λrhs  =  
      Table[ -(∂_((Subscript[y, i]))) H , 
         {i, 1, Length[λλ]}  ]

stcst  =  Flatten[ { stcst ,  λrhs } ]

Next we set up and append the rhs of system (3.3)

ξrhs  =  
      Table[ -(∂_((Subscript[z, i]))) H , 
         {i, 1, Length[ξξ]}  ]

stcst  =  Flatten[ { stcst ,  ξrhs } ]

Part Two

Vector Input For Runge Kutta Step

Vector variables and substitutions.

"st" will be the vector input to rkStep.  The order of the components of "st" is: t,y,z,λ, ξ.

st   =   
      Table[ s[ i ] ,  {i, 1 ,  1  +  n  +  m  +  n  +  m } ]

"st0" will store initial values set up later.

st0   =   
      Table[  in ,  {i, 1 ,  1  +  n  +  m  +  n  +  m } ]

In all the equations developed in Part One we must insert the correct components of "st".  To do this we make the following lists of substitutions.

ysub  =  
      Table[ yy[[i]]  ->  st[[ 1  +  i ]] , 
         {i, 1, Length[ yy ] } ]

zsub  =  
      Table[ zz[[i]]  ->  st[[ 1  + n  +  i ]] , 
         {i, 1, Length[zz]} ]

\(λsub\  = \ 
      Table[\ λλ[\([i]\)]\  -> \ 
          st[\([\ 1\ \  + n\  + m\  + \ i\ ]\)]\ , \ 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, 
          Length[\ yy\ ]}
          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

\(ξsub\  = \ 
      Table[\ ξξ[\([i]\)]\  -> \ 
          st[\([\ 1\ \  + \ n\  + m\  + n\  + \ i\ ]\)]\ , 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, 
          Length[zz]}
          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

usub  =  
      Table[ uu[[i]]  ->   nuu[[i]],  {i, 1, Length[zz]} ]

llmsub  =  
      Table[ llm[[i]]  ->   nllm[[i]]  ,  {i, 1, Length[zz]} ]

llpsub  =  
      Table[ llp[[i, i]]  ->   nllp[[i, i]]  , 
         {i, 1, Length[zz]} ]

The following, "ssub", is a master list of substitutions.

ssub  =  
      Union[ ysub ,  zsub ,  λsub ,  ξsub ,  usub , 
        llmsub ,  llpsub ]

Objects needed in vector calculations

In the following we insert the components of "st" into the objects set up in Part One.  These objects have names starting with "n".  These objects will be used in a function,  called "df", for numerically computing the rhs of the DE's during the rkStep.

nξξ  =  ξξ  //.  ssub

nF  =  F  //.  ssub

nG  =  G  //.  ssub

nG2  =  G2  //.  ssub

nDyG  =  DyG  //.  ssub

nDzG  =  DzG  //.  ssub

nstcst  =  stcst  //.  ssub

Part Three

The following module, called "df", numerically computes the right hand side of all DE's.  It uses the numerical variables set up in Part Two.  The module begins with computing the K-T multipliers and control.

\(\(\ df\ \  := \ 
      \(Module[\ {\ }\ , {

\ 
             (*\ The\ following\ solves\ for\ the\ K - 
			\ 
                T\ multiplier\ "llm".\ \ See\ paper

              				
              equations\ \((4.4)\)\ and\ \(\((4.5)\).\)\ \ \ \ \ \ \ \ \ \ *) 
              

nllm\  = \ \ 
              LinearSolve[\ \ Transpose[\ nDzG\ ]\ , \ \(-nξξ\)\ \ ]
                \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ; 


             (*\ The\ next\ line\ assigns\ the\ numerical\ values

              			
              nllm[\([i]\)]\ to\ the\ corresponding\ variable\ nlm\_i\ \ *) 

            \ Do[\ nlm\_i\  = \ nllm[\([i]\)]\ , \ {i, 1, Length[zz]}\ ]
              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ; \ 
			
\ \ \ \ \ \ \ \ 
             (*\ \ The\ following\ solves\ for\ the\ K - 
                T\ multiplier\ "llp".\ \ See\ paper\ equation\ \ 
                  \(\((4.8)\).\)\ \ *) 

            nllp\ \  = \ 
              DiagonalMatrix[\ \ 

                		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                Table[\ \ Sign[\ \ nllm[\([i]\)]\ \ ]\ \ \ nllm[\([i]\)]\ \ , 
                  
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ i, 1, Length[\ nllm\ ]\ }
                    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\ \ \ \ \ \ ]\ \ \ ; 

\ 
            Do[\ nlp\_i\  = \ nllp[\([i, i]\)]\ , \ {i, 1, Length[zz]}\ ]
              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ; 
            \ \ 
			
\ \ \ \ \ \ \ \ 
             (*\ \ The\ following\ "uurhs"\ is\ the\ right\ hand\ side\ of
                \ the\ linear\ system\ for\ the\ 

              \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              control\ "uu".\ See\ paper\ 
                \(\((4.11)\).\)\ \ \ \ \ \ \ \ \ *) 

            uurhs\  = \ 
              Table[\ Sign[\ nllm[\([i]\)]\ ] nG2[\([i]\)]\ , \ 

                  		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, Length[\ zz\ ]}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\  - 
                \ nDyG\ .\ nF\ \ \ \ ; 

\ \ \ \ \ \ \ \ \ 
             (*\ \ The\ following\ numerically\ \ solves\ \((4.11)\)\ for\ the
                \ control\ \("uu".\)\ \ *) 
\ 
            nuu\  = \ 
              LinearSolve[\ \ nDzG\ \ , \ uurhs\ \ ]
                \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ; 

\ 
            Do[\ nu\_i\  = \ nuu[\([i]\)]\ , \ {i, 1, Length[zz]}\ ]
              \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ; \ \ 

             (*\ The\ next\ instruction\ returns\ the\ numerical\ value\ of

              				the\ rhs\ of\ all\ DE' s\ \ \ \ \ *) 
					
\ 
            nstcst\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }]\)[\([1]\)]
\)\)

Part Four

Set initial values for y-state and z-state.

initialy  =  { -1.37611771672183436, 2.01517134821346122}

initialz  =  {  -1.37862693518152723, 2.01412547534158692 }

Set the initial time equal zero

st0[[1]]  =  0

Put the initial values of y and z into the initial data vector "st0".

Do[  st0[[1 + i]]  =  initialy[[i]] , 
       {i, 1, Length[ yy ]}  ]

\(Do[\ st0[\([1 + Length[\ yy\ ] + i]\)]\  = \ initialz[\([i]\)]\ , 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, Length[\ zz\ ]}\ \ ]\)

st0

Part Five

Set initial values of costate. Run Runge-Kutta solver. View graphs of output.

Set initial costate values

Adjust initial costate values until constraint set is approached.

initialλ  =  { 10.  , -10.  }

initialξ  =  { -10.  , 10. }

Put the initial values of λ and ξ into the initial data vector "st0".

\(Do[\ \ 
      st0[\([1 + Length[\ yy\ ] + Length[\ zz\ ] + i]\)]\  = \ 
        initialλ[\([i]\)]\ , 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, Length[\ yy\ ]}\ \ ]\)

\(Do[\ 
      st0[\([1 + Length[\ yy\ ] + Length[\ zz\ ] + Length[\ yy\ ] + i]\)]\  = 
        \ 
	\ \ \ \ \ \ \ initialξ[\([i]\)]\ , 
      \ \ {i, 1, Length[\ zz\ ]}
        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

st0

"tend" is the end time for solving the trajectories.  You might want to change this sometimes.

tend   =   50.

"dt" is the step size used in the Runge-Kutta step.

dt  =  .05

Run Runge-Kutta Solver

Solution is kept in table "sol".  Since this is a big table I have ended the command with a ";" so the table is not printed.  Remove the ";" to see the table.

sol  =  rk[ df,  st ,  st0,  { tend,  dt }] ;

Set up interpolating functions for graphics.

s  = 
      Flatten[ sol[[  Range[1, Length[sol] ] ,  Range[1, 1]  ]]]
         ;

\(Do[\ 
      gy[i]\  = 
        Interpolation[\ \ 
          Transpose[{\ s\ , \ 
              Flatten[\ 
						\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                sol[\([\ \ Range[1, Length[sol]\ ]\ , 

                    								\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Range[1 + i, 1 + i]\ \ ]
                    \)]]\ }]\ ]\ \ \ , 

      	\ \ \ \ \ {i, 1, Length[\ yy\ ]}
        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

\(Do[\ 
      gz[i]\  = 
        Interpolation[\ \ 
          Transpose[{\ s\ , \ 
              Flatten[\ 
						\ \ \ \ \ \ \ \ \ \ \ \ \ 
                sol[\([\ \ Range[1, Length[sol]\ ]\ \ , 

                    										\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                    Range[1 + Length[yy] + i, 
                      1 + 
										\ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
                        Length[yy] + i]\ \ ]\)]]\ }]\ ]\ \ \ , 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, Length[\ zz\ ]}
        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

\(Do[\ 
      gλ[i]\  = 
        Interpolation[\ \ 
          Transpose[{\ s\ , \ 
              Flatten[\ 
						\ \ \ \ \ \ \ \ \ \ \ \ \ 
                sol[\([\ \ Range[1, Length[sol]\ ]\ \ , 

                    										\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                    Range[1 + Length[yy] + Length[zz] + i, 
                      1 + 
										\ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
                        Length[yy] + Length[zz] + i]\ \ ]\)]]\ }]\ ]\ \ \ , 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, Length[\ yy\ ]}
        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

\(Do[\ 
      gξ[i]\  = 
        Interpolation[\ \ 
          Transpose[{\ s\ , \ 
              Flatten[\ 
						\ \ \ \ \ \ \ \ \ \ \ \ \ 
                sol[\([\ \ Range[1, Length[sol]\ ]\ \ , 

                    										\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                    Range[1 + Length[yy] + Length[zz] + Length[yy] + i, 
                      1 + 
										\ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
                        Length[yy] + Length[zz] + Length[yy] + i]\ \ ]\)]]\ }]
            \ ]\ \ \ , 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {i, 1, Length[\ zz\ ]}
        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

\(gsub\  = \ 
      Union[\ Table[\ y\_i\  -> \ \(gy[i]\)[t]\ , \ {i, 1, Length[yy]\ }\ ]
          \ , 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Table[\ z\_i\  -> \ \(gz[i]\)[t]\ , \ {i, 1, Length[zz]\ }\ ]
          \ \ \ \ \ \ ]\)

Plots of the components of the constraint G

\(\(\ Do[\ \ 
      Plot[\ Evaluate[\ G[\([i]\)]\  /. \ gsub\ ]\ , \ {\ t, 0, tend\ }\ , 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ 
          AxesOrigin\  -> \ {0, 0}\ , 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          PlotLabel\  -> \ i\ \ \ }\ ]\ \ \ \ \ \ \ \ \ , 
      \ {i, 1, Length[zz]}\ \ ]\ \ \ \ \)\)

Click the "il" to go back to initial costate conditions.

il

Compute some numerical values

tset  =  50.

\({\ Table[\ \(gy[i]\)[t]\ , \ {i, 1, Length[yy]}]\ , 
	\ 
        Table[\ \(gz[i]\)[t]\ , \ {i, 1, Length[zz]}]\ , 
	\ 
        Table[\ \(gλ[i]\)[t]\ , \ {i, 1, Length[yy]}]\ , 
	\ 
        Table[\ \(gξ[i]\)[t]\ , \ {i, 1, Length[zz]}]\ \ \ \ \ \ \ }\  /. 
      \ t\  -> \ tset\)

\(\(\ Do[\ \ 
      Plot[\ Evaluate[\ \(gλ[i]\)[t]\ ]\ , \ {\ t, 0, tend\ }\ , 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ 
          AxesOrigin\  -> \ {0, 0}\ , 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          PlotLabel\  -> \ i\ \ \ }\ ]\ \ \ \ \ \ \ \ \ , 
      \ {i, 1, Length[zz]}\ \ ]\ \ \ \ \)\)

\(\(\ Do[\ \ 
      Plot[\ Evaluate[\ \(gξ[i]\)[t]\ ]\ , \ {\ t, 0, tend\ }\ , 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\ 
          AxesOrigin\  -> \ {0, 0}\ , 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          PlotLabel\  -> \ i\ \ \ }\ ]\ \ \ \ \ \ \ \ \ , 
      \ {i, 1, Length[zz]}\ \ ]\ \ \ \ \)\)