.. _section-6.3:
Solution Techniques
-------------------
So far we have focused on the benefits derived from a generalized model,
but have not touched on the methods used to implement the model. We will
now describe the numerical techniques, first discussing the potential
problems that can occur when solving a set of equations of the form
:eq:`6.2-3` and then describing the numerical methods used to handle them.
Because the model is generalized, the solution techniques should be
transparent for a wide range of situations that can arise. In the case
of the steady-state solution finder, it is clear that the equations and
variables to be solved for are given by :eq:`6.2-4`. In certain
instances these equations may not be square yet a solution exists while
at other times the Jacobian of the right hand side may be singular. In
the case of the transient solution it is important that the solution
technique be able to maintain a user specified level of solution
accuracy under a wide range of system response times. Solution
techniques capable of handling these situations might be termed robust.
We describe such techniques here.
The reader is cautioned that the difference equations we write may
appear unconventional. The equations of reactor physics and thermal
hydraulics when differenced in one-dimensional space appear with a
single index denoting location in space. However in the generalized
control system problem, geometry or space does not enter. Instead the
continuum of space is replaced by a logical relationship among block
outputs, block inputs and plant signals. These relationships must be
recorded as part of the specific problem definition and, as will become
clear in this work, vectors can be used to store these relationships. To
summarize, the usual physical relationships that exist between points in
space in conventional difference equations are replaced by logical
relationships among blocks.
.. _section-6.3.1:
Definition of Block Diagram
~~~~~~~~~~~~~~~~~~~~~~~~~~~
The specification of the user's control equations and the solution of
these are closely tied to the block diagram concept. Therefore, before
the solution techniques can be discussed, we need to define the block
diagram.
The block diagram is represented through several vectors whose entries
define the types of blocks and the interconnections among blocks. The
vectors are one-dimensional and the index to the elements can be thought
of as analogous to a space index. To begin with, a unique signal number
is assigned to the output of a block. If there are :math:`n_{\text{blk}}`
blocks occurring in the user's input and the :math:`i`-th one is
assigned signal number m and the block is of type :math:`k` then the following
entries are created
.. math::
\begin{aligned}
L_{\text{blk,i}} = m && \text{and} && K_{\text{blk,m}} = k ; && i = 1, ..., n_{\text{blk}}~.
\end{aligned}
Further if the inputs to this block are signal :math:`q` and :math:`r` then the
entries
.. math::
\begin{aligned}
J1_{\text{blk,m}} = q && \text{and} && J2_{\text{blk,m}} = r
\end{aligned}
are also created. If the block is a non-dynamic one then the block
operator is
:math:`F_{\text{k}}\left( S_{\text{q}}^{j} , S_{\text{r}}^{j} \right)` = value of signal
:math:`m` at time :math:`j`
or, if the block is a dynamic one, then the block operator is
:math:`D_{\text{k}}\left( S_{\text{q}}^{j} , S_{\text{r}}^{j} \right)` = value of
derivative of signal :math:`m` at time :math:`j`.
The variable :math:`S_{\text{q}}^{j}` denotes that value of signal :math:`q` at time
:math:`j`. The block operators are defined in Fig. A6.1-1. An auxiliary vector
also stores information on dynamic blocks only. For the :math:`i`-th occurring
dynamic block having signal number :math:`m`, the entry
.. math::
\begin{aligned}
L_{\text{dyn,i}} = m; && i = 1, ..., n_{\text{dyn}}
\end{aligned}
is created.
A unique signal number must also be assigned to each control signal.
Recall a control signal is used to drive a plant variable and that the
signal originates at the output of a block. If there are :math:`n_{\text{ctl}}`
control signals occurring in the user's input and the :math:`i`-th one
is assigned signal number :math:`m`, then the following entries are created
:math:`L_{\text{ctl,i}} = m;\ i = 1,\ldots , n_{\text{ctl}}.`
Further, if this control signal is taken from the output of block :math:`q`,
then the entry
:math:`J_{\text{ctl,m}} = q.`
is also created.
The vectors :math:`L_{\text{blk}}`, :math:`K_{\text{blk}}`, :math:`J1_{\text{blk}}` and
:math:`J2_{\text{blk}}` thus define the block diagram. Both the
steady-state and transient solution methods access the elements in these
vectors to march through the block diagram in a manner analogous to the
step-wise progression up an axial mesh that is used in one-dimensional
thermal and hydraulic analysis codes. In the control system problem,
however, a logical relationship among blocks is substituted for the
spatial relationship among fluid cells that occurs in thermal hydraulic
problems.
.. _section-6.3.2:
Steady-State Solution
~~~~~~~~~~~~~~~~~~~~~~
The set of equations given by :eq:`6.2-4` are a set of non-linear
equations. Having the boundary conditions
:math:`{u^{*}}_{\text{mea}}\left( 0 \right)` and
:math:`{y^{*}}_{\text{ctl}}\left( 0 \right)`, we must solve this set of
equations for the initial conditions
:math:`u_{\text{dmd}}\left( 0 \right)` and :math:`x \left( 0 \right)`.
However, before describing the solution method, we discuss two important
considerations.
The first consideration is that the solution technique not require
calculation of the inverse Jacobian of :eq:`6.2-4` with respect to the
unknowns. Any non-linear equation solver that does will fail in certain
cases when in fact a solution exists. If, for example, one of the
unknowns feeds into a block that has a deadband region, then during
iteration on the unknowns the input to this block may end up in the
deadband zone. Then the derivative of the right-hand side of :eq:`6.2-4`
will be zero with respect to the unknown, in which case the inverse of
the Jacobian does not exist. The second consideration involves the
relationship between the number of equations and the number of unknowns
to be solved for. If the control system is properly designed then the
equations are either square or they are over-determined. In the later
case, the number of equations exceeds the number of unknowns. The
solution method must be able to find a solution if it exists in either
case. As for the non-linear equation solver used in this work, it will
handle both the singular Jacobian and the over-determined equation
situations. The solver is based on a least squares technique that is
described in Ref. 6-3. When used in a case with deadband as described
above, the equation solver will return a true solution if one exists.
However, the returned solution may in fact be one of an infinite number.
The non-linear equation solver uses an iterative process to converge to
the values of :math:`u_{\text{dmd}}\left( 0 \right)` and
:math:`x \left( 0 \right)` that satisfy :eq:`6.2-4`. The equation solver
provides successively more refined estimates for :math:`u_{\text{dmd}}` and
:math:`x \left( 0 \right)` at the start of every iteration and expects the
right-hand side of :eq:`6.2-4` to be calculated for each new set of
values. Values for the right-hand side are generated by marching through
the block diagram in the sequence defined by the vectors
:math:`L_{\text{blk}}`, :math:`K_{\text{blk}}`, :math:`J1_{\text{blk}}` and
:math:`J2_{\text{blk}}`. The marching procedure is as follows. Suppose
:math:`S_{\text{m}}^{P}` denotes the value of signal :math:`m` on the :math:`p`-th iterate.
Assume that if the signal :math:`S_{\text{m}}^{P}` is either a measured signal
or a control signal that it has been set with the respective boundary
condition contained in :math:`u_{\text{mea}} \left( 0 \right)` or :math:`y_{\text{ctl}} \left( 0 \right)`.
Assume also that if the signal :math:`S_{\text{m}}^{P}` corresponds to one of
the unknowns in :math:`u_{\text{dmd}} \left( 0 \right)` or :math:`x \left( 0 \right)` that the equation solver has
made its estimate for the signal prior to the :math:`p^{th}` pass
through the block diagram and the :math:`S_{\text{m}}^{P}` has been set to this
value. Then the value of the right hand-site of :eq:`6.2-4` for the
:math:`p^{th}` iteration is calculated using
.. math::
& i = 1, \ldots n_{\text{blk}} \\
& m = L_{\text{blk,i}};\ k = K_{\text{blk,m}};\ q = J1_{\text{blk,m}};\ r = J2_{\text{blk,m}} \\
& S_{\text{m}}^{P} = 0,\ k = 4 \\
& S_{\text{m}}^{P} = F_{\text{k}} \left( S_{\text{q}}^{P},\ S_{\text{r}}^{P} \right),\ k = 1,2,3,8, \ldots 21 ~.
On completion the elements of :math:`f` are
.. math::
S_{\text{J}1_{\text{blk}},\text{L}_{\text{dyn}},\text{i}}^{P}, i = 1, \ldots n_{\text{dyn}}
and the elements of :math:`g` are
.. math::
S_{\text{J}_{\text{ctl}},\text{L}_{\text{ctl}},\text{i}}^{P}, i = 1, \ldots n_{\text{ctl}} ~.
.. _section-6.3.3:
Transient Solution
~~~~~~~~~~~~~~~~~~
The numerical techniques used to solve :eq:`6.2-3` are based on explicit
differencing and a numerical marching procedure. The numerical
techniques have performed well for those problems examined to date. A
time step control mechanism automatically adjusts time step size to
maintain a specified level of accuracy. This usually results in a step
size smaller than the time constant of the fastest component which is
often a sensor. When the equations are stiff, other solution techniques
may offer a computationally more efficient solution. However, experience
has shown that the order of the control equations is usually small and
that the computational demands of the current scheme are reasonable.
The block diagram is advanced across a time step in two phases. In the
first step, the block signals are updated to the start of the time step.
This involves setting the measured and demand signals and then marching
through the block diagram while holding dynamic block signals constant
.. math::
& i = 1, \ldots n_{\text{blk}} \\
& m = L_{\text{blk,i}};\ k = K_{\text{blk,m}};\ q = J1_{\text{blk,m}};\ r = J2_{\text{blk,m}} ~, \\
& S_{\text{m}}^{j} = \frac{S_{\text{q}}^{j} - S_{\text{q}}^{j - 1}}{t^{j} - t^{j - 1}} ,\ k = 4 ~, \\
& S_{\text{m}}^{j} = F_{\text{k}}\left( S_{\text{q}}^{j}, S_{r}^{j} \right),\ k = 1,2,3,8, \ldots 21 ~, \\
& S_{\text{m}}^{j} = S_{\text{m}}^{j},\ k = 5,6,7 ~.
Then in the second step, the block signals are advanced across the time
step.
.. math::
& i = 1,....n_{\text{blk}} \\
& m = L_{\text{blk,i}};\ k = K_{\text{blk,m}};\ q = J1_{\text{blk,m}};\ r = J2_{\text{blk,m}} ~, \\
& S_{\text{m}}^{j + 1} = \frac{S_{\text{q}}^{j} - S_{\text{q}}^{j - 1}}{t^{j} - t^{j - 1}} ,\ k = 4 ~, \\
& S_{\text{m}}^{j + 1} = F_{\text{k}} \left( S_{\text{q}}^{j + 1}, S_{\text{r}}^{j + 1} \right),\ k = 1,2,3,8, \ldots 21 \\
& S_{\text{m}}^{j + 1} = D_{\text{k}} \left( \frac{S_{\text{q}}^{j} + S_{\text{q}}^{j + 1}}{2}, \frac{S_{\text{r}}^{j}
+ S_{\text{r}}^{j + 1}}{2} \right) \left( t^{j + 1} - t^{j} \right) + S_{\text{m}}^{j},\ k = 5,6,7 ~.
On completion the elements of :math:`x^{j+1}` are stored in
.. math::
S_{\text{L}_{\text{dyn},\text{i}}}^{j + 1} i = 1, \ldots n_{\text{dyn}} ~.
An accurate and stable solution to both the control equations and the
plant equations is obtained by controlling the basic time step size
known as a subinterval. The initial size of a new subinterval is
obtained by |SAS| by extrapolating rates of change in the plant
from the previous subinterval. The control equations are advanced first
over this new subinterval according to the algorithm just described. Two
time step control mechanisms can come into effect during integration of
the control equations.
The first time step mechanism attempts to limit the error in the control
equation solution that results from numerically integrating over the
subinterval. An initial estimate for this error is made after the
integration algorithm has obtained a solution at the end of the current
subinterval. The estimate is made for each element of the vector x (i.e.
dynamic blocks) by first estimating a value at the end of the current
subinterval by linearly extrapolating the change across the previous
subinterval:
.. math::
:label: 6.3-1
S_{\text{m,e}}^{j + 1} = S_{\text{m}}^{j} + \frac{S_{\text{m}}^{j} - S_{\text{m}}^{j - 1}}{t^{j} - t^{j - 1}} \left( t^{j + 1} - t^{j} \right) ~,
where
.. math::
:label: 6.3-2
m = L_{\text{dyn,i}};\ i = 1, \ldots n_{\text{dyn}} ~.
If :math:`S_{\text{m}}^{j + 1}` is the value calculated by the integration
algorithm then the error estimate is
.. math::
e_{\text{m}}^{j} = \frac{\left| S_{\text{m}}^{j + 1} - S_{\text{m,e}}^{j + 1} \right|}{\left| S_{\text{m}}^{j} \right| + \text{F5SIG}\left( m \right)}
where :math:`\text{F5SIG} \left( m \right)` is the zero crossing parameter supplied by the user as
discussed below. The solution has converged if the quantity
:math:`e_{\text{m}}^{j}` is less than the user-supplied value for the error
criterion EPSCS. If the solution has not converged, then for purposes of
control system integration only, the subinterval is bisected into two
substeps and the control equations are again advanced over the
subinterval. The error is again computed using :eq:`6.3-2` but using the
value that resulted from the previous integration in place of
:math:`S_{\text{m,e}}^{j + 1}`. If the subinterval is still not converged, it
is again bisected so now there are four substeps in the subinterval.
This process is repeated until the error between successive iterations
as defined by :eq:`6.3-2` is less than the input value for EPSCS.
The second time step mechanism limits the relative change in the control
solution over a single subinterval. Large and unrestricted changes can
lead to instability between the control system solution and the plant
solution. After the subinterval has converged as described above, the
relative change in control signals is computed via
.. math::
:label: 6.3-3
C_{\text{m}}^{j} = \frac{\left| S_{\text{m}}^{j + 1} - S_{\text{m}}^{j} \right|}{\left| S_{\text{m}}^{j} \right| + \text{F5SIG}\left( m \right)}
where
.. math::
m = L_{\text{ctl,i}};\ i = 1, \ldots n_{\text{ctl}} ~.
where :math:`\text{F5SIG} \left( m \right)` is the zero crossing parameter whose value is supplied by
the user. For the m that gives the largest value of
:math:`C_{\text{m}}^{j}` if this value :math:`C_{\text{m}}^{j}` is greater than
the user-supplied relative change criterion EPSCPL, then the subinterval
time step is cutback so that the relative change EPSCPL is just met. The
subinterval cutback size is obtained by linear interpolation so that the
new size is the value of :math:`\Delta t` that satisfies
.. math::
:label: 6.3-4
\frac{\Delta t}{\text{EPSCPL}} = \frac{t^{j + 1} - t^{j}}{C_{\text{m}}^{j}}
If the subinterval time step is cutback, then the control system
integration starts over again using the new subinterval size.
A third time step mechanism is used to limit the relative change in the
plant solution across a subinterval. This mechanism is analogous to the
second time step mechanism and is described in Ref. 6-1.