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.
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 \(n_{\text{blk}}\) blocks occurring in the user’s input and the \(i\)-th one is assigned signal number m and the block is of type \(k\) then the following entries are created
Further if the inputs to this block are signal \(q\) and \(r\) then the entries
are also created. If the block is a non-dynamic one then the block operator is
\(F_{\text{k}}\left( S_{\text{q}}^{j} , S_{\text{r}}^{j} \right)\) = value of signal \(m\) at time \(j\)
or, if the block is a dynamic one, then the block operator is
\(D_{\text{k}}\left( S_{\text{q}}^{j} , S_{\text{r}}^{j} \right)\) = value of derivative of signal \(m\) at time \(j\).
The variable \(S_{\text{q}}^{j}\) denotes that value of signal \(q\) at time \(j\). The block operators are defined in Fig. A6.1-1. An auxiliary vector also stores information on dynamic blocks only. For the \(i\)-th occurring dynamic block having signal number \(m\), the entry
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 \(n_{\text{ctl}}\) control signals occurring in the user’s input and the \(i\)-th one is assigned signal number \(m\), then the following entries are created
\(L_{\text{ctl,i}} = m;\ i = 1,\ldots , n_{\text{ctl}}.\)
Further, if this control signal is taken from the output of block \(q\), then the entry
\(J_{\text{ctl,m}} = q.\)
is also created.
The vectors \(L_{\text{blk}}\), \(K_{\text{blk}}\), \(J1_{\text{blk}}\) and \(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.
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 \({u^{*}}_{\text{mea}}\left( 0 \right)\) and \({y^{*}}_{\text{ctl}}\left( 0 \right)\), we must solve this set of equations for the initial conditions \(u_{\text{dmd}}\left( 0 \right)\) and \(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 \(u_{\text{dmd}}\left( 0 \right)\) and \(x \left( 0 \right)\) that satisfy Eq. (6.2-4). The equation solver provides successively more refined estimates for \(u_{\text{dmd}}\) and \(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 \(L_{\text{blk}}\), \(K_{\text{blk}}\), \(J1_{\text{blk}}\) and \(J2_{\text{blk}}\). The marching procedure is as follows. Suppose \(S_{\text{m}}^{P}\) denotes the value of signal \(m\) on the \(p\)-th iterate. Assume that if the signal \(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 \(u_{\text{mea}} \left( 0 \right)\) or \(y_{\text{ctl}} \left( 0 \right)\). Assume also that if the signal \(S_{\text{m}}^{P}\) corresponds to one of the unknowns in \(u_{\text{dmd}} \left( 0 \right)\) or \(x \left( 0 \right)\) that the equation solver has made its estimate for the signal prior to the \(p^{th}\) pass through the block diagram and the \(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 \(p^{th}\) iteration is calculated using
On completion the elements of \(f\) are
and the elements of \(g\) are
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
Then in the second step, the block signals are advanced across the time step.
On completion the elements of \(x^{j+1}\) are stored in
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 SAS4A/SASSYS‑1 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:
where
If \(S_{\text{m}}^{j + 1}\) is the value calculated by the integration algorithm then the error estimate is
where \(\text{F5SIG} \left( m \right)\) is the zero crossing parameter supplied by the user as discussed below. The solution has converged if the quantity \(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 \(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
where
where \(\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 \(C_{\text{m}}^{j}\) if this value \(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 \(\Delta t\) that satisfies
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.