Split and Merge Strategies for Solving Uncertain Equations Using Affine Arithmetic

The behaviour of systems is determined by various parameters. Due to several reasons like e. g. manufacturing tolerances these parameters can have some uncertainties. Corner Case and Monte Carlo simulations are well known approaches to handle uncertain systems. They sample the corners and random points of the parameter space, respectively. Both require many runs and do not guarantee the inclusion of the worst case. As alternatives, range based approaches can be used. They model parameter uncertainties as ranges. The simulation outputs are ranges which include all possible results created by the parameter uncertainties. One type of range arithmetic is the affine arithmetic, which allows to maintain linear correlations to avoid over-approximation. An equation solver based on affine arithmetic has been proposed earlier. Unlike many other range based approaches it can solve implicit non-linear equations. This is necessary for analog circuit simulation. For large uncertainties the solver suffers from convergence problems. To overcome these problems it is possible to split the parameter ranges, calculate the solutions separately and merge them again. For higher dimensional systems this leads to excessive runtimes as each parameter is split. To minimize the additional runtime several split and merge strategies are proposed and compared using two analog circuit examples.


INTRODUCTION
System uncertainties are introduced by different sources, e.g.manufacturing tolerances, ageing or indirect measurement of system parameters.To verify uncertain systems several methods can be used.Well known methods are Monte Carlo or Corner Case simulations.They choose sets of parameters and for every set a nominal simulation is performed as if the parameters were exact.The Monte Carlo method samples the parameter space randomly, whereas the Corner Case method takes samples from the corners of the parameter space.If the relation between the parameters and the system behaviour is non-monotonic the Corner Case approach misses the extreme values.The Monte Carlo method can detect them if enough samples are drawn.However, both methods do not guarantee the inclusion of the worst case and require a lot of simulation runs.As they use single points of the parameter space these methods are called point arithmetic.Alternatively, methods based on range arithmetics like interval or affine arithmetic can be used.They allow to describe a range in the parameter space and calculate the corresponding range in the solution space.Arithmetic operations are defined on ranges so that all values from an input range are mapped to an output range.This guarantees the inclusion of the extremal values mathematically.Interval arithmetic suffers from an effect called error explosion.To avoid this affine arithmetic was proposed [2].It maintains linear correlations during calculations and reduces the over-approximation in comparison to interval arithmetic.In literature some approaches exist to use range arithmetic for verifying uncertain systems.Methods for uncertain linear systems using affine arithmetic are described in [7], whereas interval methods for non-linear systems are shown in [5].
Equations of analog or mixed-signal circuits can contain arbitrary non-linear functions and can only be described implicitly: Parameter deviations can be caused by manufacturing tolerances or external influences like temperature.With decreasing structure sizes in microelectronics the influence of parameter deviations on the behaviour of analog circuits increases.An algorithm to solve implicit equations using affine arithmetic was proposed in [3,4].This method is limited in its convergence range.It only converges if the considered parameter deviations and non-linearities are small.To simulate larger circuits an extension of the algorithm was proposed in [6].It splits the parameter space, calculates the solution for the split range using affine arithmetic like before and merges the obtained solutions.The merged solution solves the original problem.In [6] all parameters are split and a merge is performed after each split.For higher dimensional systems this results in an explosion of the runtime, so a better strategy is needed.In this paper selective split strategies to split only selected parameters are explored.For transient simulations different merge strategies to postpone the merge operation are investigated as well.

PRELIMINARIES
Affine arithmetic is an enhancement of interval arithmetic.
In contrast to interval arithmetic it preserves linear correlations during the calculation so that over-approximation can be reduced.
Parameters and variables with deviations can be described as affine forms.They are denoted in the following with a hat symbol: All deviation symbols εi lie in the interval [−1 . . .1] and correspond to a certain source of uncertainty.In this way the reason for a variation of the circuit's behaviour can be tracked.The set N x contains all indexes of the deviation symbols which make up the affine form x. Vectors whose components are affine forms are denoted as x.
The affine form is symmetric to the central value x0.The sum of the absolute values of all xi denotes the maximal deviation from the central value and is called radius of the affine form: With the minimum and maximum of an affine form it can be converted to an interval: The arithmetic operations on affine forms are defined in such a way that they return an affine form.All operations on affine forms f (p) : are defined inclusion isotonely concerning the corresponding point-arithmetic function That means the following holds: All points within the input range are mapped on points of the output range.Functions can be classified in affine and non-affine functions.The result of an affine function can be exactly described as an affine form.Addition, subtraction and multiplication with a scalar constant belong to this class: All other operations are non-affine.To represent their results as an affine form a first order Taylor approximation with a Lagrangian remainder around the central value x0 is performed: The term x0 + f (x0)(x − x0) can be written as an affine form.The approximation error is enclosed by an additional deviation term yn+1εn+1 which holds the maximal approximation error in ξ ∈ [min(x), max(x)].So the result of a non-affine form can be written as: The additional deviation symbol εn+1 is uncorrelated to all previously used symbols.With each evaluation of a nonaffine function an additional deviation symbol is generated.These will be called NLPD symbols (non-linear partial deviation) in the following.The second derivatives of most elementary functions are monotone.In these cases no search for the extremal values has to be performed to find the maximal approximation error.Functions composed from several elementary functions have non-monotone deviations in general.
To avoid a computationally intensive search for the extremal values the optimal minimal inclusion is set aside.Instead, composed functions are decomposed into a sequence of elementary functions which is processed step by step.This results in an output range which overestimates the minimal range but guarantees to include it.This is one reason for over-approximation, the difference between the optimal and the calculated affine solution.Other sources of overapproximation are the solving algorithm itself and the merge operation.This will be discussed below.
An algorithm to solve implicit equations with uncertain parameters using affine arithmetic was described in [3,4].In Fig. 1 it is shown together with the extension from [6].It determines the nominal solution first and then adds deviation symbols.The deviations are calculated from the linearisation of the parameter dependencies.Due to the linearisation error the inclusion is not guaranteed.To maintain inclusion additional uncorrelated deviation symbols are added to each variable.Their size is adjusted by an iterative method.The nominal solution is calculated by the function solveNominal using the Newton-Raphson method.The affine parameters p are replaced by the nominal central values p0.The solutions x0 are used as the central values of the affine solutions x.As a first approximation of the deviation symbols xi the linearised equation system is used: These deviation symbols have the same indexes as the parameter deviations and describe the same source of uncertainty.They are called PPD symbols (parameter partial deviation).
Due to the linearisation error the solution for all ti do 8: if ti,0 > 0 then 9: si = max( ti) 10: else 11: si = − min( ti) 12: end if 13: end for 14: With the current estimation x the residuum of f ( x) is calculated.This is used to calculate a refinement of the EPD symbols.These steps are repeated until the scaling factor si for every variable xi reaches one or the xEP D,i disappears.This guarantees that the solution includes the minimal solution area under all combinations of parameter variations in p.
This algorithm can be used for all simulation types like AC, DC, transient (TR) and reachability analyses.In this paper we focus on DC and transient analysis.
The original algorithm suffers from some convergence problems.If the parameter deviations are too large or the equations are strongly non-linear the algorithm does not converge.To overcome these problems the parameter range can be split if necessary.For each part the affine solution is computed as before.Afterwards the solutions are merged to form an inclusion isotone solution of the original problem.
The split of a single parameter pi into two parts p i and p i is performed as follows: However, the split into more than two parts can be performed in a similar way.It cannot be determined a priori if it is necessary to perform a split to achieve convergence.In practice there are some critical operating points, but most operating points can be solved without splitting.Three criteria to automate the decision when to split have been proposed and compared in [6].A split is performed if the algorithm does not converge in nmax steps (Criterion 1).If the considered system is solvable the scaling factor s over the number of iterations falls monotonically and approaches 1 asymptotically.Strongly increasing values (Criterion 2) or local maxima (Criterion 3) can be additional indicators for a non-converging system.After a split the function solveAffine is called recursively with the split parameters p and p .The solutions x and x obtained from these calls are merged.The merged solution x solves the circuit's equations for both p and p as well as the original set of parameters p inclusion isotonely.The merge algorithm used in [6] is shown in Fig. 2. To obtain the results of this paper an improved version of the merge operation shown in Fig. 3 was used.It gives the same result as the original version in Fig. 2 for non-overlapping affine forms, but it avoids over-approximation if the split solutions x and x have a large overlap or are nearly identical (see Fig. 4).This happens at points with large non-linearities.Nevertheless, the inclusion of the solution is maintained.

SPLIT STRATEGIES
One possibility to overcome the limitation of the affine solving algorithm is to split the parameter ranges.The smaller the ranges are the better the algorithm converges.As we have an implicit equation system it is not possible to split in state space as proposed e. g. in [1,5].We have to split the parameter range and calculate the state space solution for that split.Splitting all parameters like proposed in [6] leads to a large computational effort for high-dimensional systems as 2 n parts have to be calculated separately. 2 n is the number for splitting n parameter in halves, the number of parts for splitting in more parts is calculated analogically.
To avoid this effort we propose to select one parameter to split.The selection is done for every split operation so different parameters can be chosen in subsequent split operations.After nmax splits using the selected split strategy were performed, the algorithm falls back to split all parameters to avoid endless loops if the current split strategy selects a parameter which does not improve convergence (see Fig. 5).
mergeSolutions( x , x ) 1: top = max(max( x ), max( x )) 2: bottom = min(min( x ), min( x )) 3: x.center = bottom+top  As the relation between parameters and the result space is non-linear in general, the optimal parameter to split can not be obtained directly.Several strategies are described and compared here.They use the helper function findMax which returns the column index of the largest element ai,j in the matrix A. The strategy "select by deviation" selects the parameter with the largest deviation: selectByDeviation 1: deviationM atrix = [p1.deviation() . . .pn.deviation()] 2: return findMax(deviationM atrix).
The strategy "select by sensitivity" selects the parameter with the largest sensitivity according to the largest entry in the sensitivity matrix df dp : selectBySensitivity 1: return findMax df dp .
Using the strategy "select by sensitivity and deviation" the sensitivity is scaled by the corresponding parameter deviations.The operator means element-wise multiplication.The parameter with the largest entry in the scaled sensitivity matrix A is selected: selectBySensitivityAndDeviation The algorithms "select by EPD and PPD" and "select by EPD and sensitivity" use the size of the EPD symbol which is adjusted during solving.Both select the equation with the largest EPD symbol.As this equation can depend on mul- return selectByX() 4: else 5: return All 6: end if Figure 5: Selection of parameter to split, select-ByX with X ∈ {Deviation, Sensitivity, Sensitivity-AndDeviation, EPDAndPPD, EPDAndSensitivity} solveAffine(solutionT ree, parameterT ree) 1: for all nodes n ∈ parameterT ree do 2: n.isSolved = f alse 3: end for 4: for all leaves l ∈ parameterT ree && l.isSolved do 5: p = l 6: replaceDerivatives() . . .

13:
l.isSolved = true 14: solutionT ree(l) = x 15: end if 16: end for In contrast, "select by EPD and sensitivity" selects the parameter with the largest sensitivity:

MERGE STRATEGIES
The algorithm in Fig. 1 always performs a merge operation after a split.For transient simulations other strategies are possible and are investigated in the following.For this purpose the solveAffine method has been modified.The pa-1: prevSolution = 0 2: prevP arameter = p 3: for all tn ∈ [0 . . .t end ] do 4: (solutionT ree, parameterT ree) = solveAffine(prevSolution, prevP arameter) 5: (prevSolution, prevP arameter) = mergeTree(solutionT ree, parameterT ree) 6: # » solution.append(prevSolution)7: end for 8: plot( # » solution) rameters and the solutions are stored in trees (see Fig. 9).Each node of the parameter tree represents a combination of parameters.The flag isSolved indicates that this combination has been solved successfully.The result for this combination is stored in the corresponding node of the solution tree.The extended solving algorithm in Fig. 6 resets the isSolved flag for all nodes (Line 1) of the parameter tree at first.Then it iterates over all leaves not yet solved (Line 4).With the parameter set stored in the current leaf the calculations according to Fig. 1 are performed.The function replaceDerivatives replaces the derivatives with the backward Euler formula.The values of the previous variables are taken from the corresponding node of the solution tree of the previous time step tn−1.If an additional split operation is necessary in the current time step tn which has not been calculated in the previous time step tn−1 the result from the corresponding parent node is taken (see Fig. 10).This introduces over-approximation as the split parameter range gives tighter inclusions in general but maintains the inclusion of the result.Otherwise it would be necessary to go back in time and recalculate all splits additionally used in the current time step for all previous time steps.If the solution step was successful the node is marked as solved and the result is stored in the solution tree (Lines 13 f.).If a split is needed then one or all parameters are selected and a split is performed.The resulting subtree is attached to the current leaf l and the iteration over all leaves of the parameter tree is restarted (Line 9 -11).After solving the equations with the split parameter a merge operation is performed to get a single affine result for plotting.
For transient simulations there are several options when the merge operation can be performed.It can be performed after each time step.This strategy is called "merge each step" (see Fig. 7).This way only time steps which require splitting according to the used splitting criteria are split.On the other hand, it is necessary to redetermine the required split parameter set for each time step.This requires lots of iterations of the affine solving algorithm as an a priori determination is not possible.As the merge operation introduces additional over-approximation, this method provides a more pessimistic inclusion.The error propagates from time step to time step via the backward Euler formula.The opposed strategy to keep the split parameter for all remaining time steps ("do not merge", Fig. 8) gives tighter inclusion as the corresponding node from the previous solution tree is used for integration.Some compromises between both extremes 1: prevSolutionT ree.root = 0 2: prevP arameterT ree.root = p 3: for all tn ∈ [0 . . .t end ] do 4: (prevSolutionT ree, prevP arameterT ree) = solveAffine(prevSolutionT ree, prevP arameterT ree) 5: # » solutionT ree.append(prevSolutionT ree) 6: end for 7: Figure 9: Parameter and solution tree built by algorithm in Fig. 6 are considered, too.The strategy "merge every n th step" keeps the splits from previous steps and only merges after n steps where n can be set by the user.Another strategy tries to solve without splits first and only if this fails it reconsiders the splits from the previous time step (strategy "try root first").

RESULTS
The affine solving algorithm and the proposed split and merge strategies were implemented in MATLAB.To compare their performances two exemplary circuits were investigated.The first one is the inverter circuit using a bipolar transistor depicted in Fig. 11.Its parameters and their uncertainties are given in Table 1.It was chosen because of its non-linear behaviour when sweeping the operating ranges of  The parameters and their uncertainties are given in Table 2.
To determine if splitting is necessary a maximum number of 30 iterations (Criterion 1) was used, unless noted otherwise.The simulation setups used for these examples are summarized in Table 3.The over-approximation and the runtimes for setup S1 with different numbers of splits (forced independently of split criterion) is shown in Fig. 13.This setup cannot be solved without splitting.The over-estimation decreases with the number of splits but the runtime increases exponentially.The number of splits and with it the runtime needed to solve the system equations increase with the size of the deviations (see Fig. 14).The time needed to solve using the same number of splits and the same size of deviations, can be reduced by using the different split strategies (see Table 4).For comparison the results of splitting a random parameter and splitting all parameters are also given.The results in Fig. 13 and Table 4 show the behaviour of the split strategies by an artificial enforcement of a given number of splits.The runtimes for the different     15.They are evaluated using DC simulations for the inverter (S1) and transient simulations for the bandpass circuit (S3) as no splitting was necessary for DC.The runtimes are normalized to the case when all parameters are split.For both circuits the strategy "select by deviation" is the fastest.The benefit is a lot larger for the bandpass (3 % of the runtime necessary) as for the inverter (74 %) as it has more uncertain parameters.For the BJT inverter the strategies "sensitivity", "sensitivity and deviation" and "random" are even slower than splitting all parameters.
For the evaluation of the merge strategies transient simulations using the setups S2 and S3 were performed.The split strategy was set to split all parameters in each split operation.As a compromise for both circuits n = 4 was selected.The runtimes for the different merge strategies are shown in Fig. 16.They are given normalized to the strategy with the largest runtime which was "do not merge" for both circuits.The fastest merge strategy for the inverter was "merge every n th time step" whereas for the bandpass "merge each step" was the fastest.A comparison of the overapproximation for the different merge strategies is shown in Fig. 17.The over-approximation was calculated as mean relative over-approximation mOA: dA and dMC denote the interval width determined from the affine and the minimum/maximum of 1000 Monte-Carlo runs, respectively.For the bandpass the over-approximation for all strategies is larger than for the BJT inverter as there are more uncertain parameters and more split operations are performed.The "do not merge" strategy causes the smallest over-approximation for both circuits.For the BJT inverter the benefit is small compared to the other three strategies which are almost equal in terms of over-approximation.By choosing the slowest strategy "do not merge" the overapproximation can be reduced most efficiently.

CONCLUSION
The convergence range of the original affine solving algorithm is limited to small deviations.If the parameters are split it is possible to perform simulations on circuits with uncertain parameters which are not possible without splitting.Additionally, the over-approximation can be reduced by employing more splits.For higher dimensional systems splitting all parameters leads to excessive runtimes as the number of parts to calculate increases exponentially.To avoid this drawback different strategies for splitting and merging were presented and compared.The split strategies select one parameter to split by different criteria.The best split strategy in terms of runtime is "select by deviation" which selects the parameter with the largest deviation for splitting.This reduces the runtime to 3 % of the time needed to split all parameters for the bandpass example.For transient simulations different merge strategies were proposed.
The fastest merge strategy depends on the circuit.For the BJT inverter the strategy "merge every n th step" and for the bandpass "merge each step" was the fastest.Comparing the merge with the split strategies it can be observed that the latter provide higher benefits in terms of runtime.Nevertheless, merge strategies are crucial as they influence the overapproximation via the backward Euler formula and by that the solvability of transient simulations.The merge strategy "do not merge" caused the smallest over-approximation but needs the largest runtime.The other merge strategies are faster but cause larger over-approximation which is almost equal among themselves.By choosing a merge strategy runtime can be traded for over-approximation depending on which effect is more critical in the current application.

Figure 6 :
Figure 6: Affine solving algorithm using trees to store results and parameters

Figure 8 :
Figure 8: Affine transient simulation with merge strategy "Do not merge"

Figure 10 :
Figure 10: Solution tree for two time steps of a transient simulation, the dotted arrows show the source node for the backward Euler formula

Figure 13 :Figure 14 :
Figure 13: Over-approximation and runtime versus number of splits for S1, split strategy = All

Figure 17 :
Figure 17: Over-approximation mOA for different merge strategies

Table 2 :
Parameters for circuit in Fig.12