Stabilization of periodic systems

PeriodicSystems.pclqr — Function

 pclqr(psys, Q, R, S; intpol = true, kwargs...) -> (F, EVALS)

Compute the optimal periodic stabilizing gain matrix F(t), such that for a continuous-time periodic state-space model psys of the form

  .
  x(t) = A(t)x(t) + B(t)u(t) + Bw(t)w(t)
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t),

the state feedback control law

 u(t) = F(t)x(t)

minimizes the quadratic index

 J = Integral {x'(t)Q(t)x(t) + u(t)'R(t)u(t) + 2*x'(t)S(t)u(t)} dt.

For a system of order n with m inputs, Q(t) and R(t) are n×n and m×m symmetric matrices, respectively, and S(t) is an n×m matrix. The matrix S is set to zero when omitted. The dimension of u(t) is deduced from the dimension of R(t). R, Q and S can be alternatively provided as constant real matrices.

Also returned are the closed-loop characteristic exponents EVALS of A(t)+B(t)F(t).

For intpol = true (default), the resulting periodic gain F(t) is computed from the stabilizing solution of a continuous-time periodic matrix differential Riccati equation using interpolation based formulas. If intpol = false, the gain F(t) is computed from a multi-point solution of the Riccati differential equation by the integration of the corresponding ODE using the nearest point value as initial condition. This option is not recommended to be used jointly with symplectic integrators, which are used by default.

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.prcric (excepting intpol).

Note: The pair (A(t),B(t)) must be stabilizable, R must be positive definite and [Q S;S' R] must be nonnegative definite .

source

PeriodicSystems.pclqry — Function

 pclqry(psys, Q, R, S; intpol = true, kwargs...) -> (F, EVALS)

Compute the optimal periodic stabilizing gain matrix F(t), such that for a periodic continuous-time state-space model psys of the form . x(t) = A(t)x(t) + B(t)u(t) + Bw(t)w(t) y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t) ,

the state feedback control law

 u(t) = F(t)x(t)

minimizes the quadratic index

 J = Integral {y'(t)Q(t)y(t) + u'(t)R(t)u(t) + 2*y'(t)S(t)u(t)} dt

For a system with m control inputs u(t) and p outputs, Q(t) and R(t) are p×p and m×m symmetric matrices, respectively, and S(t) is an p×m matrix. The dimension of u(t) is deduced from the dimension of R(t). R, Q and S can be alternatively provided as constant real matrices. The matrix S is set to zero when omitted.

Also returned are the closed-loop characteristic exponents EVALS of A(t)+B(t)F(t).

For intpol = true (default), the resulting periodic gain F(t) is computed from the stabilizing solution of a continuous-time periodic matrix differential Riccati equation using interpolation based formulas. If intpol = false, the gain F(t) is computed from a multi-point solution of the Riccati differential equation by the integration of the corresponding ODE using the nearest point value as initial condition. This option is not recommended to be used jointly with symplectic integrators, which are used by default.

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.prcric (excepting intpol).

Note: The pair (A(t),B(t)) must be stabilizable and [Q S;S' R] must be nonnegative definite.

source

PeriodicSystems.pdlqr — Function

 pdlqr(psys, Q, R, S; kwargs...) -> (F, EVALS)

Compute the optimal periodic stabilizing gain matrix F(t), such that for a discrete-time periodic state-space model psys of the form

x(t+1) = A(t)x(t) + B(t)u(t) + Bw(t)w(t) 
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t)

the state feedback control law

 u(t) = F(t)x(t)

minimizes the quadratic index

 J = Sum {x'(t)Q(t)x(t) + u'(t)R(t)u(t) + 2*x'(t)S(t)u(t)}

For a system of order n(t) with m control inputs in u(t), Q(t) and R(t) are n(t)×n(t) and m×m symmetric matrices, respectively, and S(t) is an n(t)×m matrix. The matrix S is set to zero when omitted. The dimension of u(t) is deduced from the dimension of R(t). R, Q and S can be alternatively provided as constant real matrices.

Also returned are the closed-loop characteristic exponents EVALS of A(t)+B(t)F(t).

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.prdric.

Note: The pair (A(t),B(t)) must be stabilizable and [Q S;S' R] must be nonnegative definite.

source

PeriodicSystems.pdlqry — Function

 pdlqry(psys, Q, R, S; kwargs...) -> (F, EVALS)

Compute the optimal periodic stabilizing gain matrix F(t), such that for a periodic discrete-time state-space model psys of the form

x(t+1) = A(t)x(t) + B(t)u(t) + Bw(t)w(t)
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t),

the state feedback control law

 u(t) = F(t)x(t)

minimizes the quadratic index

 J = Sum {y'(t)Q(t)y(t) + u'(t)R(t)u(t) + 2*y'(t)S(t)u(t)}

For a system with m control inputs u(t) and p outputs, Q(t) and R(t) are p×p and m×m symmetric matrices, respectively, and S(t) is a p×m matrix. The dimension of u(t) is deduced from the dimension of R(t). R, Q and S can be alternatively provided as constant real matrices. The matrix S is set to zero when omitted.

Also returned are the closed-loop characteristic exponents EVALS of A(t)+B(t)F(t).

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.prdric.

Note: The pair (A(t),B(t)) must be stabilizable and [Q S;S' R] must be nonnegative definite.

source

PeriodicSystems.pckeg — Function

 pckeg(psys, Qw, Rv, Sn; intpol = true, kwargs...) -> (L, EVALS)

Compute the Kalman estimator periodic gain matrix L(t) for a continuous-time periodic state-space model psys of the form

  .
  x(t) = A(t)x(t) + B(t)u(t) + w(t)
  y(t) = C(t)x(t) + D(t)u(t) + v(t)

and the noise covariance data E{w(t)w(t)'} = Qw(t), E{v(t)v(t)'} = Rv(t), E{w(t)v(t)'} = Sn(t), for a Kalman estimator

   .
  xe(t) = A(t)xe(t) + B(t)u(t) + L(t)(y(t)-C(t)xe(t)-D(t)u(t))

For a system of order n with p outputs, Qw and Rv are n×n and p×p symmetric matrices, respectively, and Sn is an n×p matrix. Qw, Rv and Sn can be alternatively provided as constant real matrices. The matrix Sn is set to zero when omitted.

Also returned are the closed-loop characteristic exponents EVALS of A(t)-L(t)C(t).

For intpol = true (default), the resulting periodic gain L(t) is computed from the stabilizing solution of a continuous-time periodic matrix differential Riccati equation using interpolation based formulas. If intpol = false, the gain L(t) is computed from a multi-point solution of the Riccati differential equation by the integration of the corresponding ODE using the nearest point value as initial condition. This option is not recommended to be used jointly with symplectic integrators, which are used by default.

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.pfcric (excepting intpol).

Note: The pair (A(t),C(t)) must be detectable and [Qw Sn;Sn' Rv] must be nonnegative definite.

source

PeriodicSystems.pckegw — Function

 pckegw(psys, Qw, Rv, Sn; intpol = true, kwargs...) -> (L, EVALS)

Compute the Kalman estimator periodic gain matrix L(t) for a continuous-time periodic state-space model psys of the form

  .
  x(t) = A(t)x(t) + B(t)u(t) + Bw(t)w(t)
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t) + v(t)

and the noise covariance data E{w(t)w(t)'} = Qw(t), E{v(t)v(t)'} = Rv(t), E{w(t)v(t)'} = Sn(t), for a Kalman estimator

   .
  xe(t) = A(t)xe(t) + B(t)u(t) + L(t)(y(t)-C(t)xe(t)-D(t)u(t))

For a system with mw disturbance inputs and p outputs, Qw and Rv are mw×mw and p×p symmetric matrices, respectively, and Sn is an mw×p matrix. Qw, Rv and Sn can be alternatively provided as constant real matrices. The matrix Sn is set to zero when omitted. The number of disturbance inputs mw is defined by the order of matrix Qw. Also returned are the closed-loop characteristic exponents EVALS of A(t)-L(t)C(t).

For intpol = true (default), the resulting periodic gain L(t) is computed from the stabilizing solution of a continuous-time periodic matrix differential Riccati equation using interpolation based formulas. If intpol = false, the gain L(t) is computed from a multi-point solution of the Riccati differential equation by the integration of the corresponding ODE using the nearest point value as initial condition. This option is not recommended to be used jointly with symplectic integrators, which are used by default.

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.pfcric (excepting intpol).

Note: The pair (A(t),C(t)) must be detectable, Qw must be non-negative definite, Rv must be positive definite and [Qw Sn; Sn' Rv] must be nonnegative definite.

source

PeriodicSystems.pdkeg — Function

 pdkeg(psys, Qw, Rv, Sn; kwargs...) -> (L, EVALS)

Compute the Kalman estimator periodic gain matrix L(t) for a discrete-time periodic state-space model psys of the form

x(t+1) = A(t)x(t) + B(t)u(t) + w(t)
  y(t) = C(t)x(t) + D(t)u(t) + v(t)

and the noise covariance data E{w(t)w(t)'} = Qw(t), E{v(t)v(t)'} = Rv(t), E{w(t)v(t)'} = Sn(t), for a Kalman estimator

  xe(t+1) = A(t)xe(t) + B(t)u(t) + L(t)(y(t)-C(t)xe(t)-D(t)u(t))

For a system of order n with p outputs, Qw and Rv are n×n and p×p symmetric matrices, respectively, and Sn is an n×p matrix. Qw, Rv and Sn can be alternatively provided as constant real matrices. The matrix Sn is set to zero when omitted.

Also returned are the closed-loop characteristic exponents EVALS of A(t)-L(t)C(t).

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.pfdric.

Note: The pair (A(t),C(t)) must be detectable and [Qw Sn;Sn' Rv] must be nonnegative definite.

source

PeriodicSystems.pdkegw — Function

 pdkegw(psys, Qw, Rv, Sn; kwargs...) -> (L, EVALS)

Compute the Kalman estimator periodic gain matrix L(t) for a discrete-time periodic state-space model psys of the form

x(t+1) = A(t)x(t) + B(t)u(t) + Bw(t)w(t)
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t) + v(t)

and the noise covariance data E{w(t)w(t)'} = Qw(t), E{v(t)v(t)'} = Rv(t), E{w(t)v(t)'} = Sn(t), for a Kalman estimator

  xe(t+1) = A(t)xe(t) + B(t)u(t) + L(t)(y(t)-C(t)xe(t)-D(t)u(t))

For a system with mw disturbance inputs and p outputs, Qw and Rv are mw×mw and p×p symmetric matrices, respectively, and Sn is an mw×p matrix. Qw, Rv and Sn can be alternatively provided as constant real matrices. The matrix Sn is set to zero when omitted. The number of disturbance inputs mw is defined by the order of matrix Qw.

Also returned are the closed-loop characteristic exponents EVALS of A(t)-L(t)C(t).

The keyword arguments contained in kwargs are those used for the function PeriodicMatrixEquations.pfdric.

Note: The pair (A(t),C(t)) must be detectable, Qw must be non-negative definite, Rv must be positive definite and [Qw Sn; Sn' Rv] must be nonnegative definite.

source

PeriodicSystems.pcpofstab_sw — Function

pcpofstab_sw(psys, ts = missing; K = 100, vinit, optimizer, maxit, vtol, Jtol, gtol, show_trace) -> (Fstab,info)

For a continuous-time periodic system psys = (A(t), B(t), C(t), D(t)) of period T determine a periodic output feedback gain matrix Fstab(t) of the same period and switching times ts, such that the characteristic exponents Λ of the closed-loop state-matrix A(t)+B(t)*Fstab(t)*inv(I-D(t)*Fstab(t))*C(t) are stable. The matrices of the system psys are of type PeriodicFunctionMatrix. The ns switching times contained in the vector ts must satisfy 0 = ts[1] < ts[2] < ... < ts[ns] < T. If ts = missing, then ts = [0] is used by default (i.e., constant output feedback).

The output feedback gain Fstab(t) is computed as Fstab(t) = inv(I+F(t)D(t))*F(t), with F(t) defined as

 F(t) = F_i  for t ∈ [ts[i],ts[i+1]) and i ∈ {1, ..., ns-1} or 
 F(t) = F_ns for t ∈ [ts[ns],T)

where F_i is the i-th gain. The resulting periodic matrix Fstab(t) is of type PeriodicSwitchingMatrix. The corresponding closed-loop periodic system can be obtained using the function psfeedback.

For the determination of the optimal feedback gains F_i for i = 1, ...., ns an optimization-based approach is employed using tools available in the optimization package Optim.jl. By default, the gradient-free Nelder-Mead local search method for unconstrained minimizations is employed using the keyword argument setting optimizer = Optim.NelderMead(). The alternative gradient-free Simulated Annealing global search method can be selected with optimizer = Optim.SimulatedAnnealing().

For a system with m inputs and p outputs, an internal optimization variable v is used, formed as an m×p×ns array with v[:,:,i] := F_i, for i = 1, ..., ns. The performance index to be minimized is J := sdeg(v), where sdeg(v) is the stability degree defined as the largest real part of the characteristic exponents of Af(t) := A(t)+B(t)*F(t)*C(t). The keyword argument K is the number of factors used to express the monodromy matrix of Af(t) (default: K = 100), when evaluating the characteristic exponents. By default, v is initialized as v = 0 (i.e., a zero array of appropriate dimensions). The keyword argument vinit = v0 can be used to initialize v with an arbitrary m×p×ns array v0.

The optimization process is controlled using several keyword parameters. The keyword parameter maxit can be used to specify the maximum number of iterations to be performed (default: maxit = 1000). The keyword argument vtol can be used to specify the absolute tolerance in the changes of the optimization variable v (default: vtol = 0). The keyword argument Jtol can be used to specify the relative tolerance in the changes of the optimization criterion J (default: Jtol = 0), while gtol is the method specific main convergence tolerance (default: gtol = 1e-3). With the keyword argument show_trace = true, a trace of the optimization algorithm's state is shown on stdout (default show_trace = false). (see the documentation of the Optim.jl package for additional information).

The returned named tuple info contains (fopt, sdeg0, sdeg, vopt, optres), where:

info.fopt is the resulting value of the optimal performance J;

info.sdeg0 is the initial stability degree of the closed-loop characteristic exponents;

info.sdeg is the resulting stability degree of the closed-loop characteristic exponents;

info.vopt is the resulting value of the optimization variable v;

info.optres is the result returned by the Optim.optimize(...) function of the Optim.jl package; several functions provided by this package can be used to inquire various information related to the optimization results (see the documention of this package).

source

PeriodicSystems.pcpofstab_hr — Function

pcpofstab_hr(psys,  nh = 0; K = 100, vinit, optimizer = "local", maxiter, vtol, Jtol, gtol, show_trace) -> (Fstab,info)

For a continuoous-time periodic system psys = (A(t), B(t), C(t), D(t)) of period T determine a periodic output feedback gain matrix Fstab(t) of the same period, such that the characteristic exponents Λ of the closed-loop state-matrix A(t)+B(t)*Fstab(t)*inv(I-D(t)*Fstab(t))*C(t) are stable. The matrices of the system psys are of type HarmonicArray.

The resulting output feedback gain Fstab(t) is computed as Fstab(t) = inv(I+F(t)D(t))*F(t), with F(t) in the harmonic representation form

              nh
 F(t) = F0 +  ∑ ( Fc_i*cos(i*t*2*π/T)+Fs_i*sin(i*2*π*t/T) ) ,
             i=1

where F0 is the constant term, Fc_i is the i-th cosinus coefficient matrix and Fs_i is the i-th sinus coefficient matrix. By default, the number of harmonics is nh = 0 (i.e., constant output feedback is used). The resulting periodic matrix Fstab(t) is of type HarmonicArray. The corresponding closed-loop periodic system can be obtained using the function psfeedback.

For the determination of the optimal feedback gains F0, Fc_i and Fs_i for i = 1, ...., nh an optimization-based approach is employed using using tools available in the optimization package Optim.jl. By default, the gradient-free Nelder-Mead local search method for unconstrained minimizations is employed using the keyword argument setting optimizer = Optim.NelderMead(). The alternative gradient-free Simulated Annealing global search method can be selected with optimizer = Optim.SimulatedAnnealing().

For a system with m inputs and p outputs, an internal optimization variable v is used, formed as an m*p*(2*nh+1) dimensional vector v := [vec(F0); vec(Fc_1); vec(Fs_1), ... ; vec(Fc_nh); vec(Fs_nh)]. The performance index to be minimized is J := sdeg(v), where sdeg(v) is the stability degree defined as the largest real part of the characteristic exponents of Af(t) := A(t)+B(t)*F(t)*C(t). The keyword argument K is the number of factors used to express the monodromy matrix of Af(t) (default: K = 100), when evaluating the characteristic exponents. By default, v is initialized as v = 0 (i.e., a zero array of appropriate dimensions). The keyword argument vinit = v0 can be used to initialize v with an arbitrary m*p*(2*nh+1) array v0.

The optimization process is controlled using several keyword parameters. The keyword parameter maxiter = maxit can be used to specify the maximum number of iterations to be performed (default: maxit = 1000). The keyword argument vtol can be used to specify the absolute tolerance in the changes of the optimization variable v (default: vtol = 0). The keyword argument Jtol can be used to specify the relative tolerance in the changes of the optimization criterion J (default: Jtol = 0), while gtol is the method specific main convergence tolerance (default: gtol = 1e-3). With the keyword argument show_trace = true, a trace of the optimization algorithm's state is shown on stdout (default show_trace = false). (see the documentation of the Optim.jl package for additional information).

The returned named tuple info contains (fopt, sdeg0, sdeg, vopt, optres), where:

info.fopt is the resulting value of the optimal performance J;

info.sdeg0 is the initial stability degree of the closed-loop characteristic exponents;

info.sdeg is the resulting stability degree of the closed-loop characteristic exponents;

info.vopt is the resulting value of the optimization variable v;

info.optres is the result returned by the Optim.optimize(...) function of the Optim.jl package; several functions provided by this package can be used to inquire various information related to the optimization results (see the documention of this package).

source

PeriodicSystems.pdpofstab_sw — Function

pdpofstab_sw(psys, ns = missing; vinit, optimizer, maxit, vtol, Jtol, gtol, show_trace) -> (Fstab,info)

For a discrete-time periodic system psys = (A(t), B(t), C(t), D(t)) determine a periodic output feedback gain matrix Fstab(t) of the same period, such that the characteristic exponents Λ of the closed-loop state-matrix A(t)+B(t)*Fstab(t)*inv(I-D(t)*Fstab(t))*C(t) are stable. The matrices of the system psys are of type PeriodicArray. The switching times for the resulting switching periodic gain Fstab(t) are specified by the N-dimensional integer vector ns. By default, ns = [N], where N is the maximal number of samples (i.e., N = psys.period/psys.Ts), which corresponds to a constant gain.

The output feedback gain Fstab(t) is computed as Fstab(t) = inv(I+F(t)D(t))*F(t), with F(t) defined as

F(t) = F_i for t ∈ [ns[i]Δ,ns[i+1]Δ) and i ∈ {1, ..., N-1}, or
F(t) = F_N for t ∈ [ns[N]Δ,T),

where T is the system period (i.e., T = psys.period), Δ is the system sampling time (i.e., Δ = psys.Ts), and F_i is the i-th gain. The resulting periodic matrix Fstab(t) is of type SwitchingPeriodicArray. The corresponding closed-loop periodic system can be obtained using the function psfeedback.

For the determination of the optimal feedback gains F_i for i = 1, ...., N an optimization-based approach is employed using using tools available in the optimization package Optim.jl. By default, the gradient-free Nelder-Mead local search method for unconstrained minimizations is employed using the keyword argument setting optimizer = Optim.NelderMead(). The alternative gradient-free Simulated Annealing global search method can be selected with optimizer = Optim.SimulatedAnnealing().

For a system with m inputs and p outputs, an internal optimization variable v is used, defined as an m×p×N array. By default, v is initialized as v = 0 (i.e., a zero array of appropriate dimensions). The keyword argument vinit = v0 can be used to initialize v with an arbitrary m×p×N array v0.

The performance index to be minimized is J := sdeg(v), where sdeg(v) is the stability degree defined as the largest modulus of the characteristic exponents of Af(t) := A(t)+B(t)*F(t)*C(t). The optimization process is controlled using several keyword parameters. The keyword parameter maxit can be used to specify the maximum number of iterations to be performed (default: maxit = 1000). The keyword argument vtol can be used to specify the absolute tolerance in the changes of the optimization variable v (default: vtol = 0). The keyword argument Jtol can be used to specify the relative tolerance in the changes of the optimization criterion J (default: Jtol = 0), while gtol is the method specific main convergence tolerance (default: gtol = 1e-3). With the keyword argument show_trace = true, a trace of the optimization algorithm's state is shown on stdout (default show_trace = false). (see the documentation of the Optim.jl package for additional information).

The returned named tuple info contains (fopt, sdeg0, sdeg, vopt, optres), where:

info.fopt is the resulting value of the optimal performance J;

info.sdeg0 is the initial stability degree of the closed-loop characteristic exponents;

info.sdeg is the resulting stability degree of the closed-loop characteristic exponents;

info.vopt is the resulting value of the optimization variable v;

info.optres is the result returned by the Optim.optimize(...) function of the Optim.jl package; several functions provided by this package can be used to inquire various information related to the optimization results (see the documention of this package).

source

PeriodicSystems.pdpofstab_hr — Function

pdpofstab_hr(psys, nh = 0; vinit, optimizer, maxit, vtol, Jtol, gtol, show_trace) -> (Fstab,info)

For a discrete-time periodic system psys = (A(t), B(t), C(t), D(t)) of period T determine a periodic output feedback gain matrix Fstab(t) of the same period, such that the characteristic exponents Λ of the closed-loop state-matrix A(t)+B(t)*Fstab(t)*inv(I-D(t)*Fstab(t))*C(t) are stable. The matrices of the system psys are of type HarmonicArray.

The resulting output feedback gain Fstab(t) is of type PeriodicMatrix and is computed as Fstab(t) = inv(I+F(t)D(t))*F(t), with F(t), of type PeriodicMatrix, built by sampling, with the sample time Δ = abs(psys.Ts), the harmonic representation form

              nh
Fh(t) = F0 +  ∑ ( Fc_i*cos(i*t*2*π/T)+Fs_i*sin(i*2*π*t/T) ) ,
             i=1

where F0 is the constant term, Fc_i is the i-th cosinus coefficient matrix and Fs_i is the i-th sinus coefficient matrix. F(t) is defined as F(t) = Fh((Δ(i-1)))' for t ∈ [Δ(i-1), Δi), i = 1, ..., T/Δ. By default, the number of harmonics is nh = 0 (i.e., constant output feedback is used). The corresponding closed-loop periodic system can be obtained using the function psfeedback.

For the determination of the optimal feedback gains F0, Fc_i and Fs_i for i = 1, ...., nh an optimization-based approach is employed using using tools available in the optimization package Optim.jl. By default, the gradient-free Nelder-Mead local search method for unconstrained minimizations is employed using the keyword argument setting optimizer = Optim.NelderMead(). The alternative gradient-free Simulated Annealing global search method can be selected with optimizer = Optim.SimulatedAnnealing().

For a system with m inputs and p outputs, an internal optimization variable v is used, formed as an m*p*(2*nh+1) dimensional vector v := [vec(F0); vec(Fc_1); vec(Fs_1), ... ; vec(Fc_nh); vec(Fs_nh)]'. The performance index to be minimized isJ := sdeg(v), wheresdeg(v)is the stability degree defined as the largest modulus of the characteristic exponents ofAf(t) := A(t)+B(t)F(t)C(t). By default,vis initialized asv = 0(i.e., a zero array of appropriate dimensions). The keyword argumentvinit = v0can be used to initializevwith an arbitrarymp(2*nh+1)arrayv0`.

The optimization process is controlled using several keyword parameters. The keyword parameter maxiter = maxit can be used to specify the maximum number of iterations to be performed (default: maxit = 1000). The keyword argument vtol can be used to specify the absolute tolerance in the changes of the optimization variable v (default: vtol = 0). The keyword argument Jtol can be used to specify the relative tolerance in the changes of the optimization criterion J (default: Jtol = 0), while gtol is the method specific main convergence tolerance (default: gtol = 1e-3). With the keyword argument show_trace = true, a trace of the optimization algorithm's state is shown on stdout (default show_trace = false). (see the documentation of the Optim.jl package for additional information).

The returned named tuple info contains (fopt, sdeg0, sdeg, vopt, optres), where:

info.fopt is the resulting value of the optimal performance J;

info.sdeg0 is the initial stability degree of the closed-loop characteristic exponents;

info.sdeg is the resulting stability degree of the closed-loop characteristic exponents;

info.vopt is the resulting value of the optimization variable v;

info.optres is the result returned by the Optim.optimize(...) function of the Optim.jl package; several functions provided by this package can be used to inquire various information related to the optimization results (see the documention of this package).

source

PeriodicSystems.pclqofc_sw — Function

pclqofc_sw(psys, Q, R, ts = missing; K = 1, sdeg = 0, G = I, vinit, optimizer, stabilizer,
           maxiter, vtol, Jtol, gtol, show_trace, solver, reltol, abstol, 
           N = 128, quad = false ) -> (Fopt,info)

Compute the optimal periodic stabilizing gain matrix Fopt(t), such that for a continuous-time periodic state-space model psys of the form

  .
  x(t) = A(t)x(t) + B(t)u(t) + Bw(t)w(t)  
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t) ,

the output feedback control law

u(t) = Fopt(t)*y(t),

minimizes the expectation of the quadratic index

         ∞
 J = E{ Int [x(t)'*Q(t)*x(t) + u(t)'*R(t)*u(t)]dt },
        t=0

where Q(t) and R(t) are periodic weighting matrices. The matrices of the system psys are of type PeriodicFunctionMatrix. For a system of order n with m control inputs in u(t) and p measurable outputs in y(t), Q(t) and R(t) are n×n and m×m symmetric periodic matrices of type PeriodicFunctionMatrix, respectively. The dimension m of u(t) is deduced from the dimension of R(t). Q and R can be alternatively provided as constant real matrices.

The resulting m×p periodic output feedback gain Fopt(t) is of type PeriodicSwitchingMatrix, with the switching times defined by the vector ts. The ns switching times contained in the vector ts must satisfy 0 = ts[1] < ts[2] < ... < ts[ns] < T, where T is the system period. If ts = missing, then ts = [0] is used by default (i.e., constant output feedback).

The output feedback gain Fopt(t) is computed as Fopt(t) = inv(I+F(t)D(t))*F(t), with F(t) defined as

 F(t) = F_i  for t ∈ [ts[i],ts[i+1]) and i ∈ {1, ..., ns-1} or 
 F(t) = F_ns for t ∈ [ts[ns],T)

where F_i is the i-th gain.

The covariance matrix of the initial state x(0) can be specified via the keyword argument G (default: G = I) and a desired stability degree of the closed-loop characteristic exponents can be specified using the keyword argument sdeg (default: sdeg = 0).

For the determination of the optimal feedback gains F_i for i = 1, ...., ns an optimization-based approach is employed using tools available in the optimization package Optim.jl. By default, the gradient-based limited-memory quasi-Newton method (also known as L-BFGS) for unconstrained minimizations is employed using the keyword argument setting optimizer = LBFGS(;alphaguess = LineSearches.InitialStatic(;scaled=true)), where an initial step length for the line search algorithm is chosen using the keyword argument alphaguess (see the LineSearches.jl package for alternative options). The employed default line search algorithm is HagerZhang() and an alternative method can be specified using the keyword argument linesearch (e.g., linesearch = LineSearches.MoreThuente()). Alternative gradient-based methods can be also selected, such as, for example, the quasi-Newton method BFGS with optimizer = Optim.BFGS(), or for small size optimization problems, the Nelder-Mead gradient-free method with optimizer = Optim.NelderMead(). For the computation of the function J and its gradient ∇J, the formulas developed in [1] for stable systems are used. Each evaluation involves the solution of of a pair of periodic Lyapunov differential equations using single or multiple shooting methods proposed in [2]. If the original system psys is unstable, the computation of a stabilizing feedback is performed using the same optimization techniques applied iteratively to systems with modified the state matrices of the form A(t)-αI, where α ≥ 0 is chosen such that A(t)-αI is stable, and the values of α are successively decreased until the stabilization is achieved. The optimization method for stabilization can be independently selected using the keyword argument stabilizer, with the default setting stabilizer = LBFGS(;alphaguess = LineSearches.InitialStatic(;scaled=true)). If only stabilization is desired, then use optimizer = nothing.

An internal optimization variable v is used, formed as an m×p×ns array with v[:,:,i] := F_i, for i = 1, ..., ns. By default, v is initialized as v = 0 (i.e., a zero array of appropriate dimensions). The keyword argument vinit = v0 can be used to initialize v with an arbitrary m×p×ns array v0.

The optimization process is controlled using several keyword parameters. The keyword parameter maxiter = maxit can be used to specify the maximum number of iterations to be performed (default: maxit = 1000). The keyword argument vtol can be used to specify the absolute tolerance in the changes of the optimization variable v (default: vtol = 0). The keyword argument Jtol can be used to specify the relative tolerance in the changes of the optimization criterion J (default: Jtol = 0), while gtol can be used to specify the absolute tolerance in the gradient ∇J, in infinity norm (default: gtol = 1e-5). With the keyword argument show_trace = true, a trace of the optimization algorithm's state is shown on stdout (default show_trace = false). For stabilization purposes, the values Jtol = 1.e-3, gtol = 1.e-2, maxit = 20 are used to favor faster convergence.

The returned named tuple info contains (fopt, sdeg0, sdeg, vopt, optres), where:

info.fopt is the resulting value of the optimal performance J;

info.sdeg0 is the initial stability degree of the closed-loop characteristic exponents;

info.sdeg is the resulting stability degree of the closed-loop characteristic exponents;

info.vopt is the resulting value of the optimization variable v;

info.optres is the result returned by the Optim.optimize(...) function of the Optim.jl package; several functions provided by this package can be used to inquire various information related to the optimization results (see the documention of this package).

Several keyword arguments can be used to control the integration of the involved ODE's for the solution of periodic differential Lyapunov equations for function and gradient evaluations.

If K = 1 (default), the single shooting method is employed to compute periodic generators [1]. If K > 1, the multiple-shooting method of [2] is employed, first, to convert the continuous-time periodic Lyapunov differential equations into discrete-time periodic Lyapunov equations satisfied by the generator solution in K grid points and then to compute the solution by solving an appropriate discrete-time periodic Lyapunov equation using the periodic Schur method of [3]. If quad = true, a quadrature-based evaluation of gradients is used, as proposed in [1], in conjunction with interpolation techniques. The number of sample values to be used for interpolation can be specified with the keyword parameter N (deafult: N = 128).

The ODE solver to be employed to convert the continuous-time problem into a discrete-time problem can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4), absolute accuracy abstol (default: abstol = 1.e-7) and stepsize dt (default: dt = 0, only used if solver = "symplectic"). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

Parallel computation of the matrices of the discrete-time problem can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

[1] L. Vigano, M. Bergamasco, M. Lovera, and A. Varga. Optimal periodic output feedback control: a continuous-time approach and a case study. Int. J. Control, Vol. 83, pp. 897–914, 2010.

[2] A. Varga. On solving periodic differential matrix equations with applications to periodic system norms computation. Proc. IEEE CDC/ECC, Seville, 2005.

[3] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

source

PeriodicSystems.pclqofc_hr — Function

pclqofc_hr(psys, Q, R, nh = 0; K = 1, sdeg = 0, G = I, vinit, optimizer, stabilizer,
           maxiter, vtol, Jtol, gtol, show_trace, solver, reltol, abstol, 
           N = 128, quad = false ) -> (Fopt,info)

Compute the optimal periodic stabilizing gain matrix Fopt(t), such that for a continuous-time periodic state-space model psys of the form

  .
  x(t) = A(t)x(t) + B(t)u(t) + Bw(t)w(t)  
  y(t) = C(t)x(t) + D(t)u(t) + Dw(t)w(t) ,

the output feedback control law

u(t) = Fopt(t)*y(t),

minimizes the expectation of the quadratic index

         ∞
 J = E{ Int [x(t)'*Q(t)*x(t) + u(t)'*R(t)*u(t)]dt },
        t=0

where Q(t) and R(t) are periodic weighting matrices. The matrices of the system psys are of type HarmonicArray. For a system of order n with m control inputs in u(t) and p measurable outputs in y(t), Q(t) and R(t) are n×n and m×m symmetric periodic matrices of type HarmonicArray, respectively. The dimension m of u(t) is deduced from the dimension of R(t). Q and R can be alternatively provided as constant real matrices.

The resulting m×p periodic output feedback gain Fopt(t) is of type HarmonicArray and is computed as Fopt(t) = inv(I+F(t)D(t))*F(t), with F(t) in the harmonic representation form

              nh
 F(t) = F0 +  ∑ ( Fc_i*cos(i*t*2*π/T)+Fs_i*sin(i*2*π*t/T) ) ,
             i=1

where T is the system period, F0 is the constant term, Fc_i is the i-th cosinus coefficient matrix and Fs_i is the i-th sinus coefficient matrix. By default, the number of harmonics is nh = 0 (i.e., constant output feedback is used).

The covariance matrix of the initial state x(0) can be specified via the keyword argument G (default: G = I) and a desired stability degree of the closed-loop characteristic exponents can be specified using the keyword argument sdeg (default: sdeg = 0).

For the determination of the optimal feedback gains F0, Fc_i and Fs_i for i = 1, ...., nh an optimization-based approach is employed using tools available in the optimization package Optim.jl. By default, the gradient-based limited-memory quasi-Newton method (also known as L-BFGS) for unconstrained minimizations is employed using the keyword argument setting optimizer = LBFGS(;alphaguess = LineSearches.InitialStatic(;scaled=true)), where an initial step length for the line search algorithm is chosen using the keyword argument alphaguess (see the LineSearches.jl package for alternative options). The employed default line search algorithm is HagerZhang() and an alternative method can be specified using the keyword argument linesearch (e.g., linesearch = LineSearches.MoreThuente()). Alternative gradient-based methods can be also selected, such as, for example, the quasi-Newton method BFGS with optimizer = Optim.BFGS(), or for small size optimization problems, the Nelder-Mead gradient-free method with optimizer = Optim.NelderMead(). For the computation of the function J and its gradient ∇J, the formulas developed in [1] for stable systems are used. Each evaluation involves the solution of of a pair of periodic Lyapunov differential equations using single or multiple shooting methods proposed in [2]. If the original system psys is unstable, the computation of a stabilizing feedback is performed using the same optimization techniques applied iteratively to systems with modified the state matrices of the form A(t)-αI, where α ≥ 0 is chosen such that A(t)-αI is stable, and the values of α are successively decreased until the stabilization is achieved. The optimization method for stabilization can be independently selected using the keyword argument stabilizer, with the default setting stabilizer = LBFGS(;alphaguess = LineSearches.InitialStatic(;scaled=true)). If only stabilization is desired, then use optimizer = nothing.

An internal optimization variable v is used, formed as an m*p*(2*nh+1) dimensional vector v := [vec(F0); vec(Fc_1); vec(Fs_1), ... ; vec(Kc_nh); vec(Ks_nh)]'. By default,vis initialized asv = 0(i.e., a zero vector of appropriate dimension). The keyword argumentvinit = v0can be used to initializevwith an arbitrary vectorv0`.