%% About Sector Bounds and Sector Indices
%% Conic Sectors
% In its simplest form, a conic sector is the 2-D region delimited by two lines, y=au and y=bu.
%
% The shaded region is characterized by the inequality (y-au)(y-bu)<0. More generally, any such sector can be parameterized as:
%     \left(\begin{array}{c}y\\u\end{array}\right)^T Q
% \left(\begin{array}{c}y\\u\end{array}\right) < 0,
% where Q is a 2x2 symmetric indefinite matrix (Q has one positive and one negative eigenvalue). We call Q the sector matrix. This concept generalizes to higher dimensions. In an N-dimensional space, a conic sector is a set:
%   {\bf S} = \left\{ z \in {\bf R}^N : \;\; z^T Q z < 0 \right\},
% where Q is again a symmetric indefinite matrix.
%% Sector Bounds
% Sector bounds are constraints on the behavior of a system. Gain constraints and passivity constraints are special cases of sector bounds. If for all nonzero input trajectories u(t), the output trajectory z(t) = (Hu)(t) of a linear system H(s) satisfies:
%  \int_0^T z^T(t)Q z(t) dt < 0,\quad \forall T>0 ,
% then the output trajectories of H lie in the conic sector with matrix Q. Selecting different Q matrices imposes different conditions on the system's response. For example, consider trajectories y(t) = (Gu)(t) and the following values:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}0 & -I \\ -I & 0 \end{array}\right).
% These values correspond to the sector bound:
%  \int_0^T \left(\begin{array}{c}y(t)\\u(t)\end{array}\right)^T
% \left(\begin{array}{cc}0 & -I \\ -I & 0 \end{array}\right)
%  \left(\begin{array}{c}y(t)\\u(t)\end{array}\right) dt < 0,\quad \forall T>0.
% This sector bound is equivalent to the passivity condition for G(s):
%  \int_0^T y^T(t) u(t) dt > 0,\quad \forall T>0.
% In other words, passivity is a particular sector bound on the system defined by:
%  H = \left(\begin{array}{c}G \\I \end{array}\right).
%% Frequency-Domain Condition
% Because the time-domain condition must hold for all T>0, deriving an equivalent frequency-domain bound takes a little care and is not always possible. Let the following:
% $ Q = W_1^T W_1 - W_2^T W_2
% be (any) decomposition of the indefinite matrix Q into its positive and negative parts. When W_2^T H(s) is square and minimum phase (has no unstable zeros), the time-domain condition:
%  \int_0^T (Hu)(t)^T \; Q \; (Hu)(t) dt < 0,\quad \forall T>0
% is equivalent to the frequency-domain condition:
%  H(j\omega)^H Q H(j\omega) < 0 \quad \forall \omega \in \bf{R}.
% It is therefore enough to check the sector inequality for real frequencies. Using the decomposition of Q, this is also equivalent to:
%  \left\| (W_1^T H)(W_2^T H)^{-1} \right\|_\infty  < 1 .
% Note that W_2^T H is square when Q has as many negative eigenvalues as input channels in H(s). If this condition is not met, it is no longer enough (in general) to just look at real frequencies. Note also that if W_2^T H(s) is square, then it must be minimum phase for the sector bound to hold.
% This frequency-domain characterization is the basis for sectorplot. Specifically, sectorplot plots the singular values of  (W_1^T H(j\omega))(W_2^T H(j\omega))^{-1} as a function of frequency. The sector bound is satisfied if and only if the largest singular value stays below 1. Moreover, the plot contains useful information about the frequency bands where the sector bound is satisfied or violated, and the degree to which it is satisfied or violated.
% For instance, examine the sector plot of a 2-output, 2-input system for a particular sector.

rng(4,'twister');
H = rss(3,4,2); 
Q = [-5.12   2.16  -2.04   2.17
      2.16  -1.22  -0.28  -1.11
     -2.04  -0.28  -3.35   0.00
      2.17  -1.11   0.00   0.18];
sectorplot(H,Q)

% The plot shows that the largest singular value of (W_1^T H(j\omega))(W_2^T H(j\omega))^{-1} exceeds 1 below about 0.5 rad/s and in a narrow band around 3 rad/s. Therefore, H does not satisfy the sector bound represented by Q.
%
%% Relative Sector Index
% We can extend the notion of relative passivity index to arbitrary sectors. Let H(s) be an LTI system, and let:
% $ Q = W_1^T W_1 - W_2^T W_2,    \;\;\; W_1^T W_2 = 0
% be an orthogonal decomposition of Q into its positive and negative parts, as is readily obtained from the Schur decomposition of Q. The relative sector index R, or R-index, is defined as the smallest r>0 such that for all output trajectories z(t) = (Hu)(t):
%  \int_0^T z^T(t) \; (W_1^T W_1 - r^2 W_2^T W_2) \; z(t) dt < 0,\quad \forall T>0.
% Because increasing r makes W_1^T W_1 - r^2 W_2^T W_2 more negative, the inequality is usually satisfied for r large enough. However, there are cases when it can never be satisfied, in which case the R-index is R=+\infty. Clearly, the original sector bound is satisfied if and only of R \leq 1.
% To understand the geometrical interpretation of the R-index, consider the family of cones with matrix Q(r) = W_1^T W_1 - r^2 W_2^T W_2. In 2D, the cone slant angle  \theta is related to r by
%  \tan(\theta) = r { \| W_2 \| \over \| W_1 \| }
% (see diagram below). More generally, \tan(\theta) is proportional to R. Thus, given a conic sector with matrix Q, an R-index value R<1 means that we can reduce \tan(\theta) (narrow the cone) by a factor R before some output trajectory of H leaves the conic sector. Similarly, a value R>1 means that we must increase \tan(\theta) (widen the cone) by a factor R to include all output trajectories of H. This clearly makes the R-index a relative measure of how well the response of H fits in a particular conic sector.
%
% In the diagram,
% {d_1}\frac{{\left| {W_1^Tz} \right|}}{{\left\| {{W_1}}
% \right\|}},\quad {d_2}\frac{{\left| {W_2^Tz} \right|}}{{\left\| {{W_2}}
% \right\|}},\quad R = \frac{{\left| {W_1^Tz} \right|}}{{\left| {W_2^Tz}
% \right|}},
% and
% \tan \left( \theta  \right) = \frac{{{d_1}}}{{{d_2}}} =
% R\frac{{\left\| {{W_2}} \right\|}}{{\left\| {{W_1}} \right\|}}.
% When W_2^T H(s) is square and minimum phase, the R-index can also be characterized in the frequency domain as the smallest r>0 such that:
%  H(j\omega)^H \; (W_1^T W_1 - r^2 W_2^T W_2) \; H(j\omega) < 0 \quad \forall \omega \in \bf{R} .
% Using elementary algebra, this leads to:
%  R = \max_\omega \left\| (W_1^T H(j\omega))(W_2^T H(j\omega))^{-1} \right\| .
% In other words, the R-index is the peak gain of the (stable) transfer function  \Phi(s) := (W_1^T H(s))(W_2^T H(s))^{-1}, and the singular values of \Phi(jw) can be seen as the "principal" R-indices at each frequency. This also explains why plotting the R-index vs. frequency looks like a singular value plot (see sectorplot). There is a complete analogy between relative sector index and system gain. Note, however, that this analogy only holds when  W_2^T H(s) is square and minimum phase.
%% Directional Sector Index
% Similarly, we can extend the notion of directional passivity index to arbitrary sectors. Given a conic sector with matrix Q, and a direction \delta Q, the directional sector index is the largest \tau such that for all output trajectories z(t) = (Hu)(t):
%  \int_0^T z^T(t) \; (Q + \tau \delta Q) \; z(t) dt < 0,\quad \forall T>0.
% The directional passivity index for a system G(s) corresponds to:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}0 & -I \\ -I & 0 \end{array}\right) .
% The directional sector index measures by how much we need to deform the sector in the direction \delta Q to make it fit tightly around the output trajectories of H. The sector bound is satisfied if and only if the directional index is positive.
%% Common Sectors
% There are many ways to specify sector bounds. Next we review commonly encountered expressions and give the corresponding system  H and sector matrix Q for the standard form used by getSectorIndex and sectorplot:
%  \int_0^T (Hu)(t)^T Q (Hu)(t) dt < 0,\quad \forall T>0 .
% For simplicity, these descriptions use the notation:
%  \|x \|_T = \int_0^T \|x(t)\|^2 dt,
% and omit the \forall T>0 requirement.
% Passivity
% Passivity is a sector bound with:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}0 & -I \\ -I & 0 \end{array}\right).
% Gain constraint
% The gain constraint {\left\| G \right\|_\infty } < \gamma is a sector bound with:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}I & 0 \\ 0 & -\gamma^2 I \end{array}\right).
% Ratio of distances
% Consider the "interior" constraint,
%   \| y- c u \|_T < r \| u \|_T
% where c,r are scalars and y(t)=(Gu)(t). This is a sector bound with:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}I & -cI \\ -cI & (c^2-r^2)I \end{array}\right).
% The underlying conic sector is symmetric with respect to y=cu. Similarly, the "exterior" constraint,
%   \| y- c u \|_T > r \| u \|_T
% is a sector bound with:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}-I & cI \\ cI & (r^2-c^2)I \end{array}\right).
% Double inequality
% When dealing with static nonlinearities, it is common to consider conic sectors of the form
%  a u^2 < y u < b u^2,
% where y = \phi(u) is the nonlinearity output. While this relationship is not a sector bound per se, it clearly implies:
%  a  \int_0^T u(t)^2 dt < \int_0^T y(t)u(t)dt < b \int_0^T u(t)^2 dt
% along all I/O trajectories and for all T>0. This condition in turn is equivalent to a sector bound with:
%  H(s) = \left(\begin{array}{c}\phi(.) \\ 1 \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}1 & -(a+b)/2 \\ -(a+b)/2 & ab \end{array}\right).
% Product form
% Generalized sector bounds of the form:
%  \int_0^T (y(t)-K_1 u(t))^T (y(t)-K_2 u(t)) dt < 0
% correspond to:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = \left(\begin{array}{cc}2I & -(K_2+K_1^T) \\ -(K_1+K_2^T) & K_1^T K_2 + K_2^T K_1
% \end{array}\right).
% As before, the static sector bound:
%  (y-K_1 u)^T (y-K_2 u) < 0
% implies the integral sector bound above.
% QSR dissipative
% A system y=Gu is QSR-dissipative if it satisfies:
%  \int_0^T \left(\begin{array}{c}y(t)\\u(t)\end{array}\right)^T
% \left(\begin{array}{cc}\bf{Q} & \bf{S} \\ \bf{S}^T & \bf{R} \end{array}\right)
%  \left(\begin{array}{c}y(t)\\u(t)\end{array}\right) dt > 0,\quad \forall T>0.
% This is a sector bound with:
%  H(s) = \left(\begin{array}{c}G(s)\\I \end{array}\right) , \;\;\;
%    Q = -\left(\begin{array}{cc}\bf{Q} & \bf{S} \\ \bf{S}^T & \bf{R} \end{array}\right).
% Copyright 2012 The MathWorks, Inc.