text:
# Using trajectory analysis to test and illustrate microsimulation outcomes

1. Finnish Centre for Pensions, Finland
2. University of Tampere, Finland
Research article
Cite this article as: J. Salonen, H. Tikanmäki, T. Nummi; 2019; Using trajectory analysis to test and illustrate microsimulation outcomes; International Journal of Microsimulation; 12(2); 3-17. doi: 10.34196/ijm.00198

## Abstract

We propose a new data-driven way of testing and visualizing dynamic microsimulation outcome data. The proposed statistical methodology is based on trajectory analysis (Nagin, 1999), which can be used to identify several sub-populations from a population measured longitudinally. We briefly introduce the statistical basis of trajectory analysis and discuss its use in the context of microsimulation. Finally, we report our results from the Finnish microsimulation model ELSI (Tikanmäki et al., 2014; Tikanmäki et al., 2015) to illustrate the possibilities and benefits of this technique. Trajectory analysis is available in many statistical software packages (e.g., SAS, R, Stata and Mplus). We conclude that trajectory analysis is a useful tool for investigating microsimulation outcomes.

## 1. Introduction

This paper has two goals. First, we present a statistical technique called group-based trajectory analysis and demonstrate its usefulness in a microsimulation context. We use trajectory analysis to identify unknown groups or sub-populations that can yield valuable information that is not necessarily easily accessible by other means. Second, we use three examples to show how this technique, among others, can be used to validate models and to reveal possible misspecification of the microsimulation model.

Trajectory analysis has recently gained much popularity in a number of fields including psychology, criminology (Nagin and Odgers, 2010; Nagin, 2016; van der Geest et al., 2016), sociology (Hynes and Clarkberg, 2005; Don and Mickelson, 2014), education (Kokko et al., 2008), marketing (Mani and Nandkumar, 2016) and health sciences (Nummi et al., 2014; Nummi et al., 2017a). It has increasingly been used in studies on labor market attachment (Peutere et al., 2015; Nummi et al., 2017b). To our knowledge, trajectory analysis has very seldom been used in a microsimulation context.

Dynamic microsimulation results are often reported using some ex ante classification (e.g., education, age, or labor market state). The trajectory approach offers key advantages over ex ante classification, which stem from the a priori use of an economic taxonomy. The basic advantage of trajectory analysis is that it can reveal latent patterns in longitudinal data that might otherwise remain hidden, hence complementing ex ante classifications. Furthermore, the use of a formal statistical methodology has the capacity to distinguish chance variation across individuals from real differences caused by latent sub-groups (Nagin and Odgers, 2010).

In this paper, we propose two kind of trajectory analysis implementations. First, the data can be stratified by various background factors (e.g., by labor market state) and thus the latent sub-groups are to be found within the strata. This would yield information on, for example, how common the earnings trajectories (by population state) are in the stratas investigated. Second, another way would be to classify the trajectory groups by some background factor after trajectory analysis. This would yield information on how the ex ante classifier (e.g., labor market state) is divided in trajectory groups.

Population heterogeneity is an important topic in dynamic microsimulation. Credible heterogeneity in a simulated population is a desired property of any microsimulation model. Trajectory analysis is a powerful method for demonstrating this heterogeneity. It allows for the analysis of any microsimulation outcome with interesting population heterogeneity. For instance, in the Finnish ELSI model we could analyze the individual population state, pension contributions, working lives or incomes. The technique is especially suited to cohort-based analysis, but it can also handle several cohorts simultaneously.

The approach we propose has a number of advantages. First, trajectory analysis has the potential to reveal interesting sub-groups of individuals. Second, the technique of trajectory analysis with normal distribution (like other members of the exponential family) is a well-established method, with software packages readily available for microsimulation practitioners (e.g. Jones et al., 2001; Leisch, 2004; Grun and Leisch, 2007; Haughton et al., 2009; Muthén and Muthén, 2010). Third, trajectory analysis is a flexible method that can be applied to both cross-sectional data (single period measurements), and more importantly, to longitudinal data (multiple periods of measurements). The most common metric for indexing time is age or year. Another possible metric would be time before or since a life-event.

Descriptive analysis is often inadequate for the purposes of empirical research. Trajectory analysis can provide a more rigorous examination with time-dependent covariates or risk factors affecting trajectory group membership (Nagin, 2005; Jones and Nagin, 2007). In addition to single outcome analysis, trajectory analysis also allows for the simultaneous analysis of multiple outcomes Nagin et al., 2016). Indeed, there is now a growing body of research that uses multiple trajectory analysis (e.g. Hsu, 2015; Nummi et al., 2017b). Nummi et al. (2017b) provide an example of the simultaneous modeling of employment, education, unemployment and parental leaves using a multivariate trajectory model. These abovementioned extensions also provide interesting possibilities for microsimulation data analysis.

The method has also its disadvantages. First, trajectory analysis can only be performed with discrete-time data. Most dynamic microsimulation models are defined in discrete-time, with time-point intervals of one year (Zaidi and Rake, 2001; Li and O’Donoghue, 2013; Li et al., 2014). This leads to a second disadvantage, which can be described as state vs. event microsimulation. Trajectory analysis is easiest to implement with the calendar year as the time interval, that is, using state microsimulation models.

Trajectory analysis is often applied to normally distributed data, but it is also applicable to discrete distributions such as binomial, Poisson, multinomial, etc., making this method a useful tool for the exploratory analysis of data sets. In this paper, we show three examples of trajectory analysis in microsimulation contexts. These examples cover labor market simulation topics such as wage earnings, education and pensions. Trajectory analysis of these outcomes illustrates the technique with different data distributions.

## 2. Trajectory analysis in microsimulation

### 2.1. Trajectory analysis

Trajectory analysis is, in essence, the application of finite mixture modeling to longitudinal data. It can be used for modeling the unobserved heterogeneity of individuals measured longitudinally (e.g. Nagin, 1999; Nagin, 2005). In what follows we describe the statistical background for our application of trajectory analysis.

The statistical foundation of trajectory analysis is on finite mixture modeling. Well-established in the field of statistics, the mixture modeling method concerns modeling a statistical distribution by a mixture (or weighted sum) of distributions (see Titterington et al., 1985; McLachlan and Peel, 2000). Böhning et al. (2007) give examples of the use of finite mixture modeling in various statistical applications.

Two basic modeling approaches are called growth mixture modeling and group-based trajectory modeling. They share the same analytical objective of measuring and explaining differences across population members in their developmental course. The difference between approaches lies in the way they model individual-level heterogenity in developmental trajectories.

Growth mixture models depict the average trend of outcome and individual-specific variation around the average trend with random effects using the same parameters of change (Nagin and Odgers, 2010).

Our specific aim is to identify individuals or objects in microsimulation data with the same kind of unknown developmental profiles (trajectories or sub-groups). The microsimulation population is then splitted into several sub-populations. Let $yi=(yi1,yi2,…,yiT)′$ represent the sequence of measurements on an individual i over T periods and let $fi(yi|Xi)$ denote the marginal probability distribution of $\boldsymbol{y}_{i}$ with possible time-dependent covariates $\boldsymbol{X}_{i}$. It is assumed that $f_{i}\left (\boldsymbol{y}_{i}|\boldsymbol{X}_{i}\right)$ follows a mixture of K densities

(1) ${f}_{i}\left({\mathbit{y}}_{i}|{\mathbit{X}}_{i}\right)=\sum _{k=1}^{K}{\pi }_{k}{f}_{ik}\left({\mathbit{y}}_{i}|{\mathbit{X}}_{i}\right),\sum _{k=1}^{K}{\pi }_{k}=1with{\pi }_{k}>0,$

where $πk$ is the probability of belonging to the sub-group k and $f_{ik}\left (\boldsymbol{y}_{i}|\boldsymbol{X}_{i}\right)$ is the density for the kth sub-group. Trajectory analysis can handle discrete or continuous data. The simplest choice is to use the Bernoulli distribution {0, 1} for the mixture components $f_{ik}\left (\boldsymbol{y}_{i}|\boldsymbol{X}_{i}\right)$ : but other members of the exponential family are often applied as well. It is assumed that given kth sub-group measurements are independent. For Bernoulli mixtures we can write

(2) ${f}_{ik}\left({\mathbit{y}}_{i}|{\mathbit{X}}_{i}\right)=\prod _{t=1}^{T}{p}_{itk}^{{y}_{it}}\left(1-{p}_{itk}^{1-{y}_{it}}\right),$

where the probability $p_{itk}$ is a function of covariates $\boldsymbol{X}_{i}$. For modeling the conditional distribution of $p_{itk}$ we use the logistic regression model. For the ith individual, we can then use the equation

(3) ${p}_{itk}=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left({\mathbit{x}}_{i}^{}{\beta }_{k}\right)}{1+exp\left({\mathbit{x}}_{i}^{}{\beta }_{k}\right)},$

where $\boldsymbol{x}_{i}^{}\left (t\right)$ is the tth row of $\boldsymbol{X}_{i}$ and $\boldsymbol{\beta }_{k}$ is the parameter vector for the kth sub-group. For the analysis of continuous data, one alternative is the multivariate normal distribution

(4)

where $\boldsymbol{\mu }_{ik}$ is a function of covariates $\boldsymbol{X}_{i}$ with parameters $\boldsymbol{\beta }_{k}$ and $Σik=σk2I$ a covariance matrix within kth sub-group. Thus, the measurements are assumed to be independent within sub-group k, with the variance $\sigma _{k}^{2}$. One advantage of this assumption is that it considerably simplifies the likelihood function and thus yields a computationally lighter and more stable analysis.

For modeling the trajectory mean in time t, simple linear models are usually applied, e.g. low-degree polynomials. For our three examples we used cubic polynomial model

(5) ${\mathbit{x}}^{}{\mathbit{\beta }}_{k}={\beta }_{0k}+{\beta }_{1k}t+{\beta }_{2k}{t}^{2}+{\beta }_{3k}{t}^{3}$

to model the development within the sub-group k in time (age). The mean model can also include other time-dependent covariates or time-stable covariates (risk-factors).

For parameter estimation, Maximum Likelihood (ML) estimates can be calculated by maximizing the log-likelihood $\sum \limits_{i=1}^{N}f_{i}$ over unknown parameters $βik$ and $\sigma _{k},\ k=1,\ldots ,K$ (see e.g. Nagin, 1999; Jones et al., 2001; Jones and Nagin, 2007). In most software packages, the method used for ML estimation is the EM (Expectation and Maximization) algorithm (see Dempster et al., 1977; McLachlan and Peel, 2000). The algorithm in an iterative technique involving two steps. E step finds the expected log likelihood under current parameter estimates, the subsequent M step maximizes the expected log likelihood function. These steps are iterated until the estimates converge. When applied to trajectory analysis, the E step calculates the posterior probability for sub-group membership

(6) ${w}_{ik}=\frac{{\pi }_{k}{f}_{k}\left({y}_{i}|{X}_{i},{\varnothing }_{k}\right)}{f\left({y}_{i}|{X}_{i},{\varnothing }_{k}\right)}$

under all parameter estimates $∅^k$. The estimated sub-group probabilities are then

(7) ${\stackrel{^}{\pi }}_{k}=\frac{1}{N}\sum _{i=1}^{N}{\stackrel{^}{w}}_{ik}$

Once the model parameters have been estimated, the posterior probability estimates provide a way to assign each individual to a specific sub-group or trajectory. Individuals can be assigned to the specific sub-group in which their posterior probability is highest. From the equations, we see that the assignment of individuals into groups takes into account regression parameter estimates, the number of groups, probability distribution, membership probability estimates and time span of longitudinal dataset. These are major topics in trajectory analysis and need careful consideration from microsimulation practitioner. Nagin (2005) discusses the selection of the number of groups (p. 78–87) and related statistical information criteria (p. 63–76) and the question of conditional independence (p. 26–27).

The question of groups as real entities is discussed in Nagin and Odgers (2010) and Nagin (2016). Nagin (1999) discusses the use of group membership probabilities in the calculations, as well as the links between group membership probabilities and other time-dependent covariates besides age (or time). The nature of trajectory groups in growth mixture modeling and group-based trajectory analysis is discussed in Nagin, 2005, pp. 54–56), Nagin (2016) and Nathalie et al. (2017). Don and Mickelson (2014) discuss the good practices of trajectory analysis, especially in terms of model selection.

Note that the sub-groups revealed by trajectory analysis are not fixed constructs. They are just approximations of a more complex reality. Each individual belong to a specific group with a certain probability. It would be good if we could investigate the fit of the assumed model using the identified trajectory groups. However, since the groups are based on maximum posterior probability, this would probably introduce some correlation to within-group residuals, as the actual group could also contain individuals from other groups with a small probability or weight. The correct approach from a statistical point of view, however, would be based on the investigation, using alternative covariance structures (for example, likelihood-ratio type test statistics or some other information criterion). However, this kind of testing is impossible in Nagin’s basic model. This topic has been discussed in Nagin and Odgers (2010) and Nagin (2016).

As with any statistical analysis of dynamic microsimulation outcome data, a change in the period of observation would naturally change the results of the model. The trajectory groups could also be affected to a certain extent. However, we believe that in our examples the basic main grouping structure would remain quite similar.

### 2.2. Finnish ELSI microsimulation model

ELSI is a longitudinal microsimulation model (Tikanmäki et al., 2014; Dekkers and Van den Bosch, 2016) that is used to assess the development of the statutory pensions in Finland. The dynamic ageing model has been developed at the Finnish Centre for Pensions, a statutory co-operation body providing research and expertise services related to the Finnish pension system.

ELSI model has been designed to assess the future earnings-related pensions and the national pensions (Tikanmäki et al., 2017). It can also be used to analyze changes in the pension system and in the underlying demographic or macroeconomic conditions. One of its uses has been to assess the distributional effects of the pension reform of 2017 (Tikanmäki et al., 2015).

The Finnish pension system is mainly based on pension rights accrued on the basis of the individual’s life-time earnings. The only exception is the survivor’s pension, which accounts for no more than 6% of total pension expenditure. The ELSI model is therefore based on individual-level information and calculations of pensions received in one’s own right. The model comprises both pension recipients and those still working. The model simulates each individual’s working life prior to retirement.

The base population consists of all adults aged 18 or over covered by the social insurance system in Finland. Most of the material is drawn from administrative records maintained by the Finnish Centre for Pensions and the Social Insurance Institution of Finland. Information on educational level from Statistics Finland is also added to the model. The ELSI model provides the opportunity to run large populations. For the present analysis we used a random sample of 32% of the base population, which in numerical terms translates into 1.5 million individuals in the baseline year of 2012. The population is then simulated until 2085. Deceased people remain in the model and a new cohort as well as new immigrants enter the model each year. Consequently the population increases over the course of the simulation.

The ELSI model has a modular structure (Figure 1). Each colored box represents a module of the model, while white boxes represent external sources of information. There are no feedback loops from later to earlier modules. The simulation starts from the population module, followed by the earnings module, pension modules, and the taxation module, and finally brings together the results.

Figure 1

The population module has several functions. It simulates population and labor market transitions as well as educational changes. The population module is based on transition probabilities that are estimated from historical data for 2010–2014. The module also uses Statistics Finland’s official population forecast to replicate general trends in the sample population. Transition probabilities, with one-year time steps, are by and large deterministic, that is, based on exogenous information.

In the population module, we simulate a new population or labor market state for each individual based on transition probabilities. There are 22 states in the model, the most important of which are: active (employed), inactive, unemployed, various pension states for different pension benefits, and deceased.

The labor market transitions are of Markovian type, which means that the transition probabilities are based on current state rather than former history. However, it is possible to add memory to a Markov process by extending the state space. For instance, in model ELSI we have three different active states. One for those employed first consecutive year, another for those employed second year and third one for the rest. Hence, unemployment risks may be higher for those who do not yet have an established labor market position.

Education is also simulated in the population module. Post-basic education dynamics is based on age and gender-specific transition probabilities. Changes in education level are possible at any age, although they are not very common after age 35.

The transition probabilities are updated each simulation year using the population level information produced by the semi-aggregated LTP model (Tikanmäki et al., 2017). There is thus a simple alignment of microsimulation outcomes to macro level aggregates.

Wage earnings are simulated in the earnings module. Wages are simulated for each individual annually based on the labor market state and an underlying earnings-equation, which is a time-series model with a stochastic component. The earnings equation takes also into account gender, age and level of education. Wages are simulated for active workers and those in partial retirement. The earnings module is described in Tikanmäki et al. (2014).

The earnings-related pension module calculates pension amounts based on the simulated time of retirement (retirement age) and life-course earnings. The pension calculation takes no shortcuts but is as detailed as possible given the data in use. The national pension and guarantee pension amounts are calculated as a residual of the earnings-related pension, since they are income tested. The calculation of national pensions is described in Sihvonen (2015).

The taxation module finalizes the substantive simulation. The previous modules have produced gross wages and pensions. In the taxation module, the current (2016) rules for income taxation are applied to both simulated wages and pension earnings. The calculation of income taxes is described in Sihvonen (2015).

After the simulation run, the results are analyzed in the results module, which calculates aggregate results over the course of the simulation, based on individual-level outcomes. Measures of the distribution (mean, percentage points, Gini coefficient) of pensions can be produced by ex ante classifiers such as gender, level of education and year of birth. Aggregate measures on the duration of working life and partition of life-course into active and passive stages are also calculated in the result module. The module collects individual-level output data containing information on labor market state, wage earnings, residence, education level, pension earnings, pension benefit, working life, pension accrual, etc. Therefore the material is also available for the statistical analysis illustrated in this paper. The proposed statistical technique could be used with many other outcomes as well.

In the following analysis we use a longitudinal individual-level data set that covers the period from 2008–2085 with a 25% sub-sample of the original output data. The analysis is based on three outcomes: wage earnings, education level and pension earnings. The wage and education trajectories are presented for the male cohort born in 1995 and pension trajectories for the male cohort born in 1960. With trajectory analysis, it would be possible to analyze several cohorts at the same time so long as the data is longitudinal.

Similar longitudinal data is produced in many other dynamic microsimulation models, too. For example, the Swedish SESIM model and the Norwegian MOSART model are in this respect similar to ELSI, providing rich individual or household level labor market output data (Fredriksen and Stølen, 2007; Klevmarken and Lindgren, 2008; Flood et al., 2012).

### 2.3. Trajectory model selection

Choosing the number of mixture components K is an important stage of finite mixture modeling. The trajectory method requires the researcher to specify the assumed number of sub-groups in the data. The optimal number of sub-groups K can be assessed using information criteria, which are widely used in this context.[1] One commonly used criterion to assess the model fit is the Bayesian information criterion (BIC). The model with largest BIC is preferred (SAS implementation). One may also calculate the group-specific average posterior probability over individuals to measure the fit. If this values exceeds the minimum threshold (at least 0.7) the fit is considered satisfactory. Nagin (2005, p. 63–77) provides an overview of model selection in trajectory analysis context. Ultimately, meaningful real-world interpretations of sub-groups requires not just information criteria, but also theoretical evaluation and good judgment.

In trajectory analysis missing data is considered missing completely at random according to the taxonomy layed by Little and Rubin (1987). Software packages handle missing data in such a way that all available data is used in the estimation. In other words, all individuals or objects with some missing longitudinal data values are included in the analysis.

From Table 1 we can see the BIC values of our example analyses on three outcome-variables (earnings, education and pension). The BIC score indicates that, a five-group solution fits the wage earnings and pension variables, while for education BIC increases with more groups. However, increasing the number of education groups further would yield some infinitesimal groups. Therefore our choice is the five-group solution for the wage earnings and pension outcomes, and the six-group solution for the education outcome variable.

Table 1

In this study the computations were carried out using SAS software package with accompanying PROC TRAJ application, an easy to use PC SAS procedure for analyzing Nagin’s model (e.g. Jones et al., 2001; Jones and Nagin, 2007; Andruff et al., 2009).[2] Appendix A.1 gives an example of the programming code for wage earnings trajectories (Example 1, k=5).

The trajectory plots (Figures 24) present conditional means of time points calculated over the simulation period. Relative sizes $π^k$ of the sub-groups are also presented in the figures. These plots are the main tool for interpreting the results obtained. The estimated model is summarized in Appendix (A.2), which includes group-specific parameter estimates. The SAS procedure also produces a confidence interval (95%) and model-predicted values by sub-group if necessary.

Figure 2
Figure 3
Figure 4

It must be emphasized that the trajectories in the following examples are group means, which indicates that there is a range of individual tracks around the trajectory curves. Software packages routinely output the group estimates and confidence intervals to visualize within-group population heterogeneity. For the sake of readability, the model fit curves and confidence intervals are not presented in this paper.

In the examples we use cubic age model (see Equation 5). The choice depends on the assumed development of the data over time. One advantage in PROC TRAJ is that the age or time model can be group-specific, which means that for example one group can be fitted with a second order polynomial model, whereas another group can be fitted with a cubic model etc.

## 3. Examples

### 3.1. Earnings trajectories: binary outcome

For this example, we have chosen a young cohort to illustrate a simulated life-course. The cohort, born in 1995, enters labor markets at age 17 (year 2012), and subsequent life-course including individual wage earnings is simulated in the ELSI model’s population and earnings modules.

The example shows a trajectory analysis with binary outcome data. Individual yearly wage is dichotomized (yes = 1/no = 0), so wage earnings could also be interpreted as employment. In terms of content, wage earnings could be analyzed without transformations using continuous data (euros).

The BIC score (Table 1) indicates that a five-group solution yields the best fit for the earnings outcome, but we also present the solutions for three and four-groups solutions. Figure 2 shows that the relative sizes $π^k$ vary somewhat with an increasing number of groups. Nevertheless all solutions have essentially the same major sub-groups.

In terms of content, we can see that there is a large group (k = 3: 49%, k = 4: 33% and k = 5: 29%) with strong labor market integration, or a strong earnings profile, until retirement age. Other groups show an early declining trend in employment. Labor market integration is weak in one group (k = 3: 12%, k = 4: 10% and k = 5: 8%), especially at older age. Unemployment and permanent disability draw individuals out of employment prior to old-age retirement. Retirement age for old-age pension in this cohort is 67 years and 9 months. Partial old-age pension can be drawn three years before retirement age.

#### 3.1.1. ELSI model recalibration

This trajectory analysis experiment revealed a slight misspecification in the ELSI model’s earnings module (see spikes at age 65 in all trajectory solutions). The spikes are visible in mean-based trajectories, but they would not be seen in fitted means models. This goes to show how useful group-based mean curves can be in revealing model misspecification.

In the earnings module the mechanism for working after retirement did not work as intended. During the course of simulation, many of those who had stopped working earlier suddenly decided to return to work after drawing their pension (Figure 2). Of course, such behavior is not plausible on a larger scale.

The Finnish pension system allows full-time employment even after retirement on an old-age pension. However, wage earnings for old-age pensioners are typically quite low. The misspecification therefore did not have a major impact on the main results and was not observed in comparison with the LTP macro model.

In the ELSI model earnings are calibrated with the LTP model. In practice, this is done by comparing projected average earnings each year by gender, population state and age. Any deviations imply changes in the ELSI model’s parameters.

In the LTP model, working after retirement is not modeled in detail, and the relevant point of reference is provided by the corresponding statistical figures. In the ELSI model, the share of new retirees working after retirement is the same as in the observed statistics. Trajectory analysis showed that the pool of retirees working after retirement also included individuals who had not worked in the year preceding retirement, which was not intended.

Following the trajectory analysis of simulated earnings, the ELSI model was recalibrated appropriately. There was no change in the aggregate figures, namely wage sum and number of people working after retirement.

Trajectory analysis is not yet routinely part of ELSI model testing, but it could be incorporated as part of its validation procedures.

##### 3.1.1.1. Composition of trajectory groups

Checking the composition of the groups is a good way to validate also the trajectory analysis results. It is good practice to cross-tabulate the trajectory groups by the available background factors of the microsimulation. We have added an example which shows the labor market state or population state of earnings groups (k = 5) at age 60 in the simulated life-course (Figure 3). Depending on the available information of the microsimulation model, also other background factor could be cross-tabulated.

The labor market states (Table 2) confirm the earnings trajectory results. The high and stable wage earnings (employment) group G4 (29%) consists of people who are mainly (89.9%) in active labor market states. The weak attachment group G1 (8%) consists of individuals who had a disability or inactive status at the time. Also the share of deceased (15.5%) is higher in this group than in the other groups, which is quite natural. Another weak attachment group is G3 (17%), where the share of disabled and inactive individuals is relatively high compared to the other groups.

Table 2

### 3.2. Education trajectories: class outcome

In this second example we continue to study the young cohort born in 1995 and illustrate a simulated life-course in terms of education level. This cohort finishes compulsory education at age 16 (in 2011). The subsequent education dynamics is simulated in the education module of the ELSI model.

The example shows a trajectory analysis with class-level outcome data. The distribution of the levels of education skewed to left (low education), therefore the outcome is modeled using the zero-inflated Poisson model. Yearly individual education information records the highest level of education. In practice, each individual level of education either remains unchanged or rises during the course of the simulation. Individuals enter the model at age 18 when 99% are at the basic education level. The ELSI model includes four levels of education: 0 = basic education (ISCED 0–3), 1 = secondary education (ISCED 4), 2 = lower academic degree (ISCED 5, 5B) and 3 = higher academic degree (ISCED 5A, 6).

The BIC score (Table 1) indicates that a six-group solution yields the best fit for the education outcome, but we also present the four-group and five-group solutions. Figure 4 shows that with the six-group solution, the smallest group comprises no more than 2% of the population.

The trajectory means clearly show how education dynamics play out in practice. For example, education increases sharply until around age 30, and then is virtually stagnant at age 40. In the four-group solution, the results are dominated by the final education level. The four-group and five-group solutions show adult and further education in action, with 6–7% of the cohort still on an upward education trajectory at later ages (k = 5: 6%, k = 6: 7%).

Overall, we have observed that 20% have a low education, 40% complete secondary education, 21% complete a lower university degree and 15% a higher university degree.

### 3.3. Pension trajectories: continuous outcome

The cohort in this third example was born in 1960. For this cohort, earnings, labor market state, accrued pensions and disability pension data are covered up to age 52 (in 2012). Thereafter, pension stock and new pensions are simulated in the ELSI model’s earnings-related pension module.

The example shows a trajectory analysis with continuous normal data. The outcome variable is earnings-related pension (log euros/month), excluding survivor’s pensions, expressed in nominal value. If an individual is not retired (not yet retired or deceased) the outcome value is zero, otherwise the individual is retired and receives a pension (Figure 5). The level of pension is determined by current pension rules for different types of pension benefit, such as disability pensions, partial old-age pensions and old-age pensions.

Figure 5

Pension levels are affected by two factors. First, the most important factor that affects all retirees is indexation: earnings-related pensions are index-linked to prices (80%) and wages (20%). Second, some individuals continue to work while drawing their pension, which accrues additional pension benefits.

The BIC score (Table 1) indicates that the five-group solution yields the best fit for the pension outcome, but we also present the solutions for three-group and four-group solutions.

Trajectory analysis reveals several substantive issues. First, pensions increase as the cohort ages and approaches minimum retirement age at 64 years and 5 months (persons born in 1960 can retire flexibly between ages 64 years 5 months and 69 years). Pensions taken earlier are for the most part either disability pensions or partial old-age pensions. Second, the trajectories are affected by mortality as all trajectory means go to zero with the death of individuals in the group. For instance, one group (k = 5: 14%) dies early and its pension level remains low. This group, along with (k = 5: 11%), enters disability pensions at a higher rate than other groups. In addition, some groups have a higher life expectancy. The group with the highest life expectancy (k = 5: 18%) also has the highest pensions. Third, there is some indication of segregation among pensioners here.

It is also worth noting that pension indexation rules increase pensions in those sub-groups where life expectancy is high. Finally, an increase in the number of groups leads to the emergence of some interesting sub-groups. The early transition to disability pension is evident in the five-group solution (k = 5: 11%).

## 4. Conclusion

We have applied finite mixture modeling to longitudinal data in the context of dynamic microsimulation. It is useful in visualizing dynamic microsimulation results because it demonstrates, within a formal statistical structure, population heterogeneity in simulated datasets. Furthermore, the sub-groups identified in the analysis may reveal model misspecifications if such issues are present.

Our first goal was to illustrate the results of dynamic microsimulation, which are often presented as statistical moments or other distribution measures for a given population. Statistical moments presented on the basis of ex ante classification rules (e.g. gender, level of education, occupation or region) are a common way of reporting results for specific population sub-groups. We have shown how trajectory analysis can be applied as a data-driven method for illustrating longitudinal results in a novel way. These two ways of presenting results are mutually complementary.

Our second goal was to test the individual-level or object-level outcomes of dynamic microsimulation. Testing should be understood here by locating groups with unwanted or otherwise peculiar outcomes. Microsimulation easily yields biases and model misspecification. Microsimulation practitioners apply various techniques to calibrate models and to locate misspecification. Trajectory analysis could also be a useful tool in testing simulation model assumptions, as model misspecification may be revealed in developmental trajectories. Trajectory analysis together with other testing methods can improve microsimulation model reliability, especially in connection with large-scale revisions of model parameters.

There are some limitations with the proposed technique. First, it takes some effort to analyze all the cohorts included in the microsimulation model. It might be advisable to focus on some key cohorts if the simulation model assumes similar parameters across cohorts. Second, given the computational requirements of the EM algorithm, it is apparent that extremely large data sets (hundreds of thousands of individuals) should be avoided. A fraction or sample of the simulated data set usually reveals the underlying sub-populations with sufficient reliability. Third, it is necessary to mention conditional independence, even though this is not a serious drawback in a microsimulation context. Finally, basically trajectory analysis cannot handle nonlinear models. The group-specific regression model needs to be linear in parameters. However, reasonable approximations can often be provided by low-degree polynomial models.

Trajectory analysis is also a highly flexible method. In our case we performed the analysis for one outcome at the time, but the technique and software is in place to model several outcomes at the same time. This is called multi-trajectory modeling, which identifies groups of individuals or objects that follow similar trajectories across multiple outcomes. For example, it would be possible to analyze income, education and labor market status simultaneously. This would yield a slightly different probability for belonging to a specific group. Another possibility would be to add time-dependent (e.g. earnings) or time-stable (e.g. gender, year of birth or nationality) covariates alongside age or time polynomial. This would be a case of multinomial logistic regression within a population sub-group, and the result would consist of risk factors influencing group membership.

## Footnotes

### 1.

SAS package PROC TRAJ counts BIC and AIC routinely.There are other information criteria, like cross-validation available, if the number of sub-groups is critical for the microsimulation practitioner (see Nielsen et al., 2014).

### 2.

SAS and STATA (see Jones and Nagin, 2013) trajectory analysis package can be downloaded from the website: https://www.andrew.cmu.edu/user/bjones/. R package for trajectory analysis is Flexmix (see Leisch, 2004).

## Appendix

### A.1 An example SAS PROC TRAJ-code.

Earnings trajectories: Five-group solution

• Males born 1995

• Simulated yearly wage earnings 2015–2080

• Binary outcome 1 = some wage, 0 = no wage

• WMALE_1995 is the individual level data set with yearly observations on binary wage (bwage) and age (age).

The PROC TRAJ produces four output datasets. OP includes the group mean estimates and confidence intervals (95%). OS includes the mixture components (sizes of the sub-groups). OF includes individual-level information on group posterior probabilities and group assignment. This is an important data set where we can merge any other relevant background information. Finally, OE includes the regression model coefficients and information criteria (AIC and BIC). ITDETAIL statement shows the progression of maximum likelihood estimation in the SAS log window.

Algorithm 1

### PROC TRAJ.

1: PROC TRAJ DATA=WMALE_1995 OUTPLOT=OP OUTSTAT=OS OUT=OF OUTEST=OE ITDETAIL;

2: ID idno;

3: VAR bwage2015-bwage2080;

4: INDEP age2015-age2080;

5: MODEL logit;

6: NGROUPS 5;

7: ORDER 3 3 3 3 3;

8: RUN;

The actual model is specified with VAR, INDEP, MODEL and ORDER statements. The MODEL logit states that we analyze binary outcomes. Other possibilities are censored normal outcome (cnorm) and zero-inflated Poisson outcome (zip). NGROUP states the number of sub-groups. ID states the number of individuals. ODRER defines the number of mixture components (=NGROUPS) and respective degree of age-polynomial. Each sub-group can have different degrees of age-polynomial.

Algorithm 2

### TRAJPLOT-macro.

1: %TRAJPLOT(OP,OS,,);

2: RUN;

3: DATA WMALE_1995; MERGE WMALE_1995 OF (KEEP=idno group); BY idno;

4: RUN;

The TRAJPLOT macro statement invokes the result plot, which by default plots the trajectory means and their estimates. Finally, individual level results are collected and merged with microsimulation output dataset.

### A.2 Summary of Estimated Model: Outcomes, Groups, Parameter Estimates and their Standard Errors.

Outcome Groups $β^0$ $SE(β^0)$ $β^1$ $SE(β^1)$ $β^2$ $SE(β^2)$ $β^3$ $SE(β^3)$
Earnings 1 −0.60816 0.93616 0.17156 0.07229 −0.00845 0.00173 0.00008 0.00001
Earnings 2 10.77241 0.7512 −1.05243 0.05715 0.03048 0.00135 −0.00026 0.00001
Earnings 4 −24.0041 0.69707 1.88334 0.05421 −0.04084 0.00128 0.00026 0.00001
Earnings 5 2.89272 0.93197 −0.55508 0.08212 0.02763 0.00225 −0.00029 0.00002
Education 1 150.7634 13.69866 −23.324 1.03504 0.60943 0.01352 −0.00442 0.00018
Education 2 −38.2479 1.84202 1.95442 0.10906 −0.03343 0.00214 0.00019 0.00001
Education 3 −1.40435 0.08235 0.09431 0.00638 −0.00206 0.00016 0.00001 0.00000
Education 4 −3.82746 0.09413 0.27969 0.00712 −0.00568 0.00017 0.00004 0.00000
Education 5 −6.14041 0.09384 0.46159 0.00701 −0.00958 0.00017 0.00006 0.00000
Education 6 −34.1529 1.02762 2.02103 0.06446 −0.03841 0.00133 0.00024 0.00001
Pension 1 1352.226 1.29912 −66.1983 0.03537 1.07053 0.0006 −0.00572 0.00001
Pension 2 −220.279 3.93005 5.68868 0.15254 −0.03308 0.00196 −0.00003 0.00001
Pension 3 −464.879 3.77328 14.13735 0.14801 −0.12224 0.00197 0.00024 0.00001
Pension 4 27.60321 4.88433 −1.86771 0.19721 0.03987 0.00262 −0.00025 0.00001
Pension 5 −127.459 2.91581 3.18956 0.10992 −0.01691 0.00136 −0.00002 0.00001

## References

1. 1
Latent class growth modelling: a tutorial (2009)
Tutorials in Quantitative Methods for Psychology 5:11–24.
https://doi.org/10.20982/tqmp.05.1.p011
2. 2
Computational Statistics & Data Analysis 51:5205–5210.
https://doi.org/10.1016/j.csda.2006.10.025
3. 3
Applications of microsimulation modelling
(2016)
Prospective microsimulation of pensions in European Member States, Applications of microsimulation modelling.
4. 4
Maximum likelihood estimation for incomplete data via the EM algorithm
(1977)
Journal of the Royal Statistical Society B 39:1–38.
5. 5
Relationship satisfaction trajectories across the transition to parenthood among low-risk parents (2014)
Journal of Marriage and Family 76:677–692.
https://doi.org/10.1111/jomf.12111
6. 6
SESIM III - A Swedish dynamic micro simulation model
(2012)
SESIM III - A Swedish dynamic micro simulation model.
7. 7
Model 1: MOSART (Dynamic Cross-Sectional Microsimulation Model)
(2007)
In: Anil Gupta, Ann Harding, editors. Modelling Our Future: Population Ageing, Health and Aged Care (International Symposia in Economic Theory and Econometrics, Volume 16). Emerald Group Publishing Limited. pp. 433–437.
8. 8
FlexMix: an R package for finite mixture modelling
(2007)
Faculty of Commerce – Papers.
9. 9
Review of three latent class cluster analysis packages: latent GOLD, poLCA, and MCLUST (2009)
The American Statistician 63:81–91.
https://doi.org/10.1198/tast.2009.0016
10. 10
Trajectories of multimorbidity and impacts on successful aging (2015)
Experimental Gerontology 66:32–38.
https://doi.org/10.1016/j.exger.2015.04.005
11. 11
Women’s employment patterns during early parenthood: a group-based trajectory analysis (2005)
Journal of Marriage and Family 67:222–239.
https://doi.org/10.1111/j.0022-2445.2005.00017.x
12. 12
Advances in group-based trajectory modeling and an SAS procedure for estimating them (2007)
Sociological Methods & Research 35:542–571.
https://doi.org/10.1177/0049124106292364
13. 13
A note on a STATA plugin for estimating group-based trajectory models (2013)
Sociological Methods & Research 42:608–613.
https://doi.org/10.1177/0049124113503141
14. 14
A SAS procedure based on mixture models for estimating developmental trajectories (2001)
Sociological Methods & Research 29:374–393.
https://doi.org/10.1177/0049124101029003005
15. 15
Simulating an Ageing Population: A microsimulation approach applied to Sweden
(editors) (2008)
Bingley: Emerald Group Publishing Ltd.
16. 16
Trajectories based on Postcomprehensive and higher education: their correlates and antecedents (2008)
Journal of Social Issues 64:59–76.
https://doi.org/10.1111/j.1540-4560.2008.00548.x
17. 17
FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in R (2004)
Journal of Statistical Software 11:1–18.
https://doi.org/10.18637/jss.v011.i08
18. 18
A survey of dynamic microsimulation models: uses, model structure and methodology
(2013)
International Journal of Microsimulation 6:3–55.
19. 19
Handbook of Microsimulation Modelling. Contributions to economic analysis
(2014)
Dynamic models, Handbook of Microsimulation Modelling. Contributions to economic analysis, Emerald.
20. 20
Statistical Analysis with Missing Data
(1987)
New York: John Wiley & Sons.
21. 21
22. 22
Finite Mixture Models
(2000)
New York: John Wiley and Sons Inc.
23. 23
MPlus Statistical Analysis With Latent Variables, User’s Guide
(2010)
MPlus Statistical Analysis With Latent Variables, User’s Guide, Los Angeles, CA.
24. 24
Analyzing developmental trajectories: a semiparametric, group-based approach (1999)
Psychological Methods 4:139–157.
https://doi.org/10.1037/1082-989X.4.2.139
25. 25
Group-Based Modeling of Development
(2005)
Cambridge, MA: Harvard University Press.
26. 26
Group-Based trajectory modeling and criminal career research (2016)
Journal of Research in Crime and Delinquency 53:356–371.
https://doi.org/10.1177/0022427815611710
27. 27
The Wiley Handbook of Developmental Psychopathology
(2017)
Developmental Trajectories of Psychopathology. An Overview of Approaches and Applications, The Wiley Handbook of Developmental Psychopathology, John Wiley & Sons Ltd.
28. 28
Group-Based criminal trajectory analysis using cross-validation criteria (2014)
Communications in Statistics - Theory and Methods 43:4337–4356.
https://doi.org/10.1080/03610926.2012.719986
29. 29
Group-Based multi-trajectory modelling
(2016)
Statistical Methods in Medical Research pp. 1–9.
30. 30
Group-Based trajectory modeling (nearly) two decades later (2010)
Journal of Quantitative Criminology 26:445–453.
https://doi.org/10.1007/s10940-010-9113-7
31. 31
A trajectory analysis of body mass index for Finnish children (2014)
Journal of Applied Statistics 41:1422–1435.
https://doi.org/10.1080/02664763.2013.872232
32. 32
Trajectories of a set of ten functional somatic symptoms from adolescence to middle age (2017a)
Archives of Public Health 75:11.
https://doi.org/10.1186/s13690-017-0178-8
33. 33
New Advances in Statistics and Data Science
(2017)
Statistical Analysis of Labor Market Integration: A Mixture Regression Approach, New Advances in Statistics and Data Science, ICSA Book Series in Statistics, Springer.
34. 34
Job Contract at Birth of the First Child as a Predictor of Women’s Labor Market Attachment: Trajectory Analyses over 11 Years (2015)
Nordic Journal of Working Life Studies 5:9–30.
https://doi.org/10.19154/njwls.v5i1.4763
35. 35
Eläketurvakeskuksen ELSI-mikrosimulointimallin laajennus Kelan eläkkeisiin ja verotukseen
(2015)
Eläketurvakeskuksen ELSI-mikrosimulointimallin laajennus Kelan eläkkeisiin ja verotukseen, Finnish Centre for Pensions Working papers 03/2015.
36. 36
Microsimulating Finnish earnings-related pensions
(2014)
Finnish Centre for Pensions Working Papers 02/2014.
37. 37
Distributional effects of forthcoming Finnish pension reform – a dynamic Microsimulation approach
(2015)
International Journal of Microsimulation 8:75–98.
38. 38
Statutory pensions in Finland – long-term projections 2016
(2017)
Statutory pensions in Finland – long-term projections 2016, Finnish Centre for Pensions Reports 02/2017.
39. 39
Statistical analysis of finite mixture distributions
(1985)
Great Britain: John Wiley & Sons.
40. 40
The effects of incarceration on longitudinal trajectories of employment: a follow-up in high-risk youth from ages 23 to 32
(2016)
Crime & Delinquency 62:107–140.
41. 41
Dynamic Microsimulation Models: A Review and Some Lessons for SAGE
(2001)
The London School of Economics.

## Article and author information

### Author details

1. #### Janne Salonen

Finnish Centre for Pensions, Helsinki, Finland
##### For correspondence
janne.salonen@etk.fi
##### Competing interests
No competing interests reported
2. #### Heikki Tikanmäki

Finnish Centre for Pensions, Helsinki, Finland
##### For correspondence
heikki.tikanmaki@etk.fi
##### Competing interests
No competing interests reported
3. #### Tapio Nummi

University of Tampere, Tampere, Finland
##### Competing interests
No competing interests reported

### Acknowledgements

The authors wish to thank Statistics Finland for providing material for this study and two anonymous referees as well as chief editor Matteo Richiardi for their many thoughtful comments that greatly improved this paper.

### Publication history

1. Version of Record published: August 31, 2019 (version 1)


text:
# 1. Introduction¶

## 1.1 Problem Statement¶

This project draws inspiration from a specific aspect of a popular application called SplitWise. SplitWise basically helps people settle expenses among different people in a group in an efficient manner. It's easiest to understand with an example, so here's a minimal example (which we will actually get working later) :

• Bob pays \$20 for gas, split between him and Alice. (Alice owes Bob \$10)
• Charles pays \$20 for lunch, split between him and Bob. (Bob owes Charles \$10)

How can they settle their debts? Of course, the simplest option is :

• Alice pays Bob the \$10 she owes him • Bob pays Charles the \$10 he owes him.

If Alice, Bob and Charles are friends (we'll say they're in the same group from here on) then of course, it's easy for us to understand that intuitively, the most effective way to settle their debts would just be :

• Alice pays Charles the \$10 that Bob owes him. This reduces the total number of transanctions required for the group to reach a settled state from 2 to 1. The utility of such an application is obvious. However, what's interesting is that this problem of finding the minimum number of transanctions that will lead to a group reaching a "settled" state, is actually NP-Complete. Also, it can actually be modelled in some interesting ways, including as an integer programming problem! Over the course of this notebook, we will see one possible formulations of this problem, and also explore some extensions of the problem statement. ## 1.2 Problem Complexity¶ This seemingly straightforward problem can actually be reduced to the well known subset-sum problem, which is NP-Complete. The subset-sum problem is defined as follows : Given a set of Integers$S$, determine the existence of a non-empty subset who's elements add up to a given integer$n$One possible transformation of our problem into the subset-sum problem is as follows : Initialize the subset$S$with all amounts of money owed Set each integer$n$to the amount of money owed to each individual This means that finding an optimal solution to this problem for large sets of transanctions is computationally intractable, and interestingly enough, according to the founder of the actual SplitWise app they use a greedy solution (which, as greedy algorithms often go, isn't always optimal but run in reasonable time). ## 1.3 Data¶ The data used in this notebook is synthetically generated by a Python script that writes to a CSV file. The CSV file is in the format : (owed by), (owed to), (amount owed) This format was chosen because even though real transanctions most likely occur in the format : (person who paid), (amount paid), (people splitting the expense) Transforming from that format to the one used by this notebook is simple, and therefore does not affect the modelling in any way. It would be simple to add some pre-processing to convert data between these two formats. # 2. Mathematical Model¶ To understand the integer programming problem, let's go back to our minimal example and try and solve that. We can represent the transanctions in terms of a matrix$T$such that : each element$T_{ij}$of the matrix is the amount of money (in dollars) person$i$owes person$j$. Here's a quick recap of the transanctions we had, simplified down : Alice owes Bob \$10 Bob owes Charles \$10 So based on our above formulation of the$T$matrix:$T_{alice, bob} = 10T_{bob, charles} = 10$For convenience, let's assign integer indices to Alice, Bob and Charles as follows : Alice = 1 Bob = 2 Charles = 3 We can now form our$T$matrix :$T = \begin{bmatrix} 0 & 10 & 0 \\ 0 & 0 & 10 \\ 0 & 0 & 0 \end{bmatrix}$Note here that$T_{i,i}$will always be zero (as you can't owe yourself any money). Similarly, we can construct a settlement matrix$S$which represents the optimal solution. Let's define each element of$S$,$S_{i,j}$as they money that person$i$needs to pay person$j$in order for all debts to be settled. In our simple case, we can construct this by observation as :$S = \begin{bmatrix} 0 & 0 & 10 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$So our settlement matrix tells us that Alice must pay Charles \$10, which is what we want. This is in some sense, our "target" matrix.

Intepreted in a different way, our $S$ matrix is also an alternate to our $T$ matrix, because it's just telling us who owes who, and how much, which is exactly what the $T$ matrix too. In other words, it basically says "Alice finally owes charles \$10." The important thing to understand here, is that we don't really care who pays who, as long as everyone receives or pays what they're supposed to. We can express this as : Net Amount for person$i$= Amount owed to$i$- Amount owed by$i$How do we represent that using the matrices we've constructed? Net amount owed by$i$=$\sum_{j} T_{ij}$(sum along row$i$) Net amount owed to$i$=$\sum_{j} T_{ji}$(sum along column$i$) Since the net amount owed to or by each person must be the same in the the$S$and$T$matrices we arrive at this important constraint :$\sum_{j} T_{ji} - \sum_{j} T_{ij} = \sum_{j} S_{ji} - \sum_{j} S_{ij}$Therefore, as long as this constraint holds,$S$and$T$are basically equivalent in the manner we require. Our task has now reduced to finding a sparisified version of$T$that satisifies the above constraint! How do we actually sparisfy$T$though? One simple way of doing this is to express this as our objective: Minimize the number of elements of$S$that are not 0. And to do this, we can turn to integer programming, by simply using a binary indicator variable to determine if an element of$Sis zero or nonzero. Therefore, our final formulation is as follows : \begin{aligned} \underset{I \in \mathbb{R^{n x n}}}{\text{minimize}}\qquad& \sum_{i, j} I_{ij} && i,j = 1,\dots,n \\ \text{subject to:}\qquad& S_{i,j} \ge 0 && i,j = 1,\dots,n\\ & S_{i,j} \le MI_{i,j} && i,j = 1,\dots,n\\ & \sum_{j} S_{ji} - \sum_{j} S_{ij} = \sum_{j} T_{ji} - \sum_{j} T_{ij} && i = 1,\dots,n\\ \end{aligned} whereM$is an upper-bound on the values$S_{i,j}$can take. We can now solve an MWE of our problem. # 3. Solution¶ ## 3.1 Minimal Working Solution¶ In [1]: using JuMP, NamedArrays, Cbc  Let's reuse our example from above. In [2]: T = [0 10 0; 0 0 10; 0 0 0 ]  Out[2]: 3×3 Array{Int64,2}: 0 10 0 0 0 10 0 0 0 In [3]: m = Model(solver = CbcSolver()) @variable(m, s[1:3, 1:3]) @variable(m, indicator[1:3, 1:3], Bin) # while some possible formulations of this problem might allow # for a negative S(i,j) value in the sense that a negative # amount of money owed can be construed as money lent, but the # formulation described above does not allow for negative values of S. @constraint(m, s .>= 0) @expression(m, rowsums[i in 1:3], sum(s[i,j] for j in 1:3)) @expression(m, colsums[i in 1:3], sum(s[j,i] for j in 1:3)) # This is the constraint that forces equivalence between # the S matrix and the T matrix. @constraint(m, rowsumc[i in 1:3], (colsums[i] - rowsums[i]) == (sum(T[j,i] for j in 1:3) - sum(T[i,j] for j in 1:3))) # This constraint forces the indicator variables to be set to 1 # if a transanction occurs i.e. if any element of the S matrix is non-zero # we're just going to use an arbitrarily large number for the upper-bound for now # more on this later. @constraint(m, s .<= indicator*10000) @objective(m, Min, sum(indicator)) status = solve(m)  Out[3]: :Optimal Let's examine the value of our S matrix : In [4]: s_opt = getvalue(s)  Out[4]: 3×3 Array{Float64,2}: 0.0 0.0 10.0 0.0 0.0 0.0 0.0 0.0 0.0 This matches what we were expecting, which is great. Now that we have a simple working model, we can extend it to more realistic examples, with several transanctions. ## 3.2 Implementation¶ Now that we have a minimum working example, let's test it out with a larger data set. ### 3.2.1 Base Implementation¶ We will use a random set of 20 transanctions among 5 people. First, let's just parse the CSV file. In [5]: function get_T_from_csv(datafile) transanctions = [] all_users = [] f = open(datafile) lines = readlines(f) for line in lines owed_by, owed_to, amount = split(line, ",") push!(transanctions, [owed_by, owed_to, parse(Int, amount)]) push!(all_users, owed_by) push!(all_users, owed_to) end # Set up some variables we'll need for later unique_users = sort(unique(all_users)) n = length(unique_users) # Initialize the T Matrix T = NamedArray(zeros(n,n), (unique_users, unique_users), ("OWED BY", "OWED TO")); # Fill in the T Matrix # T[i,j] => i owes j T[i,j] for transanction in transanctions owed_by, owed_to, amount = transanction T[owed_by, owed_to] = amount end return T, n, unique_users end  Out[5]: get_T_from_csv (generic function with 1 method) In [6]: T, n, unique_users = get_T_from_csv("dense.csv") println(T)  5×5 Named Array{Float64,2} OWED BY ╲ OWED TO │ "Albert" "Anthony" … "Tyler" "Victoria" ──────────────────┼────────────────────────────────────────────────────── "Albert" │ 0.0 435.0 … 435.0 180.0 "Anthony" │ 65.0 0.0 185.0 185.0 "Elizabeth" │ 170.0 185.0 180.0 65.0 "Tyler" │ 185.0 290.0 0.0 15.0 "Victoria" │ 185.0 280.0 … 185.0 0.0  Here, we have a total of 20 transanctions that occured. Let's examine them. In [7]: for i in unique_users for j in unique_users if T[i,j] != 0 println("$i owes $j \$$(T[i,j])") end end end  Albert owes Anthony 435.0 Albert owes Elizabeth 65.0 Albert owes Tyler 435.0 Albert owes Victoria 180.0 Anthony owes Albert 65.0 Anthony owes Elizabeth 210.0 Anthony owes Tyler 185.0 Anthony owes Victoria 185.0 Elizabeth owes Albert 170.0 Elizabeth owes Anthony 185.0 Elizabeth owes Tyler 180.0 Elizabeth owes Victoria 65.0 Tyler owes Albert 185.0 Tyler owes Anthony 290.0 Tyler owes Elizabeth 290.0 Tyler owes Victoria 15.0 Victoria owes Albert 185.0 Victoria owes Anthony 280.0 Victoria owes Elizabeth 290.0 Victoria owes Tyler 185.0  Here's a graph showing us just how messy these debts are right now : Let's see how many transanctions we can reduce this down to. In [8]: function find_S_matrix(T, n, unique_users) m = Model(solver = CbcSolver()) @variable(m, s[1:n, 1:n]) @variable(m, indicator[1:n, 1:n], Bin) @constraint(m, s .>= 0) @expression(m, rowsums[i in 1:n], sum(s[i,j] for j in 1:n)) @expression(m, colsums[i in 1:n], sum(s[j,i] for j in 1:n)) @constraint(m, rowsumc[i in 1:n], (colsums[i] - rowsums[i]) == (sum(T[j,i] for j in 1:n) - sum(T[i,j] for j in 1:n))) @constraint(m, s .<= indicator*10000) @objective(m, Min, sum(indicator)) status = solve(m) println(status) s_opt = getvalue(s) S = NamedArray(s_opt, (unique_users, unique_users), ("OWED BY", "OWED TO")); return S end  Out[8]: find_S_matrix (generic function with 1 method) In [10]: S = find_S_matrix(T, n, unique_users)  Optimal  Out[10]: 5×5 Named Array{Float64,2} OWED BY ╲ OWED TO │ "Albert" "Anthony" … "Tyler" "Victoria" ──────────────────┼────────────────────────────────────────────────────── "Albert" │ 0.0 255.0 … 0.0 0.0 "Anthony" │ 0.0 0.0 0.0 0.0 "Elizabeth" │ 0.0 0.0 0.0 0.0 "Tyler" │ 0.0 0.0 0.0 0.0 "Victoria" │ 0.0 290.0 … 205.0 0.0 In [12]: for i in unique_users for j in unique_users if S[i,j] != 0 println("i owes j \$$(S[i,j])") #println("\"$i,$j,$(S[i,j])\",")
end
end
end

Albert owes Anthony $255.0 Albert owes Elizabeth$255.0
Victoria owes Anthony $290.0 Victoria owes Tyler$205.0

And here's our new, simplified graph :

Just 4 transanctions! Let's verify that this result ensures everyone is getting what they deserve.

First, let's see what everyone's receiving (we know that this is the sum along columns) :

In [13]:
colsums = []
for i in 1:n
colsum = sum(S[j,i] for j in 1:n)
# we'll also save the colsum for later
push!(colsums, colsum)
println( "$(unique_users[i]) =>$colsum" )
end

Albert => 0.0
Anthony => 545.0
Elizabeth => 255.0
Tyler => 205.0
Victoria => 0.0

Next, let's see what everyone's paying (sum along the rows):

In [14]:
rowsums = []
for i in 1:n
rowsum = sum(S[i,j] for j in 1:n)
push!(rowsums, rowsum)
println( "$(unique_users[i]) =>$rowsum" )
end

Albert => 510.0
Anthony => 0.0
Elizabeth => 0.0
Tyler => 0.0
Victoria => 495.0

Finally, let's see the net amount everyone receives or gets. This is just what they're getting minus what they owe.

In [15]:
net = colsums - rowsums
for i in 1:n
println( "$(unique_users[i]) =>$(net[i])" )
end

Albert => -510.0
Anthony => 545.0
Elizabeth => 255.0
Tyler => 205.0
Victoria => -495.0

No surprises there. Let's repeat this process for the $T$ matrix we started off with. If we end up at the same thing, we know for sure that we're being fair to everyone.

In [16]:
colsums = []
for i in 1:n
colsum = sum(T[j,i] for j in 1:n)
push!(colsums, colsum)
end

rowsums = []
for i in 1:n
rowsum = sum(T[i,j] for j in 1:n)
push!(rowsums, rowsum)
end

net = colsums - rowsums
for i in 1:n
println( "$(unique_users[i]) =>$(net[i])" )
end

Albert => -510.0
Anthony => 545.0
Elizabeth => 255.0
Tyler => 205.0
Victoria => -495.0

Great, it matches our values from the $S$ matrix.

So now we know that our model works, as no one is paying more than they owe, or receiving less than they should.

### 3.2.2 Circular Transanctions¶

An interesting case to try out would be to have many circular transanctions.

Basically, something along the lines of :

Alice owes Bob \$10 Bob owes Charles \$10
Charles owes Dick \$10 Dick owes Alice \$10

We can of course tell that in this case, there's no transanctions needed as the group is already in a settled state. But let's see if our model realizes that.

Let's get our data first

In [17]:
T, n, unique_users = get_T_from_csv("circular.csv")

Out[17]:
(10×10 Named Array{Float64,2}
OWED BY ╲ OWED TO │     "Adam"    "Amanda"  …      "Noah"     "Roger"
──────────────────┼──────────────────────────────────────────────────
"Adam"            │        0.0         0.0  …       100.0         0.0
"Amanda"          │        0.0         0.0            0.0         0.0
"Carolyn"         │        0.0         0.0            0.0         0.0
"Cheryl"          │        0.0         0.0            0.0         0.0
"Jonathan"        │        0.0         0.0            0.0         0.0
"Joshua"          │        0.0         0.0            0.0       100.0
"Juan"            │      100.0         0.0            0.0         0.0
"Nancy"           │        0.0         0.0            0.0         0.0
"Noah"            │        0.0       100.0            0.0         0.0
"Roger"           │        0.0         0.0  …         0.0         0.0, 10, SubString{String}["Adam", "Amanda", "Carolyn", "Cheryl", "Jonathan", "Joshua", "Juan", "Nancy", "Noah", "Roger"])
In [18]:
for i in unique_users
for j in unique_users
if T[i,j] != 0
println("$i owes$j \$$(T[i,j])") end end end  Adam owes Noah 100.0 Amanda owes Carolyn 100.0 Carolyn owes Nancy 100.0 Cheryl owes Joshua 100.0 Jonathan owes Juan 100.0 Joshua owes Roger 100.0 Juan owes Adam 100.0 Nancy owes Cheryl 100.0 Noah owes Amanda 100.0 Roger owes Jonathan 100.0  If we were to graph these, it becomes clear why I keep calling them "circular transanctions" : Let's give them to our solver and see what it comes up with. In [19]: S = find_S_matrix(T, n, unique_users);  Optimal  In [20]: S.array  Out[20]: 10×10 Array{Float64,2}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 It works as we expected it to. Now we can get more complex. ### 3.2.3 Multiple Groups¶ #### 3.2.3.1 Motivation¶ In a more realistic situation, there isn't going to be just one group of people on the app. There are going to exist multiple groups. Let's define what a group is exactly : • A user can belong to multiple groups. • An input transanction will only occur among members of one group. • A settlement transanction can be across groups only if it's between two people who belong to both groups. • The contrapositive to the above is that if money is flowing between two people, they must belong to at least one group together. Here's a simple example to illustrate why this is important : Example 1 Group 1 (Alice, Bob and Charles) Alice owes Bob \10 Bob owes Charles \10 Group 2 (Alice, Dick, and Charles) Charles owes Dick \20 Dick owes Alice \20 If we just treated them as individual groups, we'd end up with these 2 transanctions to stabilize each group : Alice owes Charles \10 Charles owes Alice \20 Whereas what we really want is : Charles owes Alice \10 So how do we approach this problem? We definitely can't combine everything into one big T matrix. Consider the following example to understand why : Example 2 Group 1 (Alice, Bob and Charles) Alice owes Bob \10 Bob owes Charles \10 Group 2 (Dick and Charles) Charles owes Dick \10 Here, we wouldn't want to end up with an answer that looks like : Alice owes Dick \10 Because Alice and Dick aren't in the same group together. This time, we do want our solver to give us two separate transanctions that respect the group structuring. Or in other words : Alice owes Charles \10 Charles owes Dick \10 #### 3.2.3.2 Implementation¶ It turns out, all we have to do is tell our solver who money can't move between. This is just an additional (set of) constraints that we need to derive based on the group information. But first off, we need this group information, which we will generate another CSV file for, in the format : Group Number, User1, User2... For simplicity, we will assume the group numbers start at 1 and increase by 1 with each line. Let's parse our group information CSV file. In [21]: function get_group_info(groupfile) groups = [] all_users = [] f = open(groupfile) lines = readlines(f) for line in lines this_group = split(line, ",") users_only = this_group[2:end] push!(groups, users_only) for user in users_only push!(all_users, user) end end return groups, sort(unique(all_users)) end groups, unique_users = get_group_info("groups.csv");  Here's our groups : In [22]: group_number = 1 for group in groups println("Group group_number : ") for name in group print("name, ") end println("\n") group_number+=1 end  Group 1 : Alice, Bob, Charles, Evan, Group 2 : Charles, Dick, Ron, Alice, Group 3 : Evan, Dick, Ingrid,  Now, we need to construct a list of disallowed transanctions. This is basically a list of (user1, user2) pairs that money cannot flow between. Basically, if two users are never in the same group, then money can't flow between them. In [23]: function get_disallowed_transanctions(groups, unique_users) disallowed = [] n = length(unique_users) for i in 1:n-1 for j in i:n first_user = unique_users[i] second_user = unique_users[j] allowed = false for k in 1:3 is_first_in_group_k = first_user in groups[k] is_second_in_group_k = second_user in groups[k] are_both_in_group_k = is_first_in_group_k & is_second_in_group_k if are_both_in_group_k allowed=true # no need to check the other groups as # we found at least one group in common they share break end end if allowed == false # money can't flow in either direction # for disallowed users push!(disallowed, (i,j)) push!(disallowed, (j,i)) end end end return disallowed end disallowed_transanctions = get_disallowed_transanctions(groups, unique_users);  Now that we have that, let's load up their transanctions as usual. The trasnanction file has the same structure as before. In [24]: T, n, unique_users = get_T_from_csv("group_transanctions.csv");  We need to change our model to include the additional constraint that now, money can't flow between certain people. In [27]: function find_S_matrix_grouped(T, n, unique_users) m = Model(solver = CbcSolver()) # first, all the same variables and constraints as before. @variable(m, s[1:n, 1:n]) @variable(m, indicator[1:n, 1:n], Bin) @constraint(m, s .>= 0) @expression(m, rowsums[i in 1:n], sum(s[i,j] for j in 1:n)) @expression(m, colsums[i in 1:n], sum(s[j,i] for j in 1:n)) @constraint(m, rowsumc[i in 1:n], (colsums[i] - rowsums[i]) == (sum(T[j,i] for j in 1:n) - sum(T[i,j] for j in 1:n))) @constraint(m, s .<= indicator*10000) # and now, the additional constraint of disallowed transanctions for user_tuple in disallowed_transanctions # we can actually define a constraint either on indicator, or S itself. i = user_tuple[1] j = user_tuple[2] @constraint(m, s[i,j] == 0) end # objective doesn't change either @objective(m, Min, sum(indicator)) status = solve(m) println(status) s_opt = getvalue(s) S = NamedArray(s_opt, (unique_users, unique_users), ("OWED BY", "OWED TO")); return S end;  Let's run this updated model. In [28]: S = find_S_matrix_grouped(T, n, unique_users);  Optimal  In [29]: for i in unique_users for j in unique_users if S[i,j] != 0 println("i owes j \$$(S[i,j])")
end
end
end

Charles owes Alice $350.0 Charles owes Ron$215.0
Evan owes Alice $205.0 Evan owes Bob$80.0
Evan owes Dick $230.0 Ingrid owes Dick$300.0

If we examine the results, we can actually see that none of the transanctions occur between two people who are never in the same group at least once. This is exactly what we wanted. Here's the groups again just so you can look for yourself :

In [30]:
group_number = 1
for group in groups
println("Group $group_number : ") for name in group print("$name, ")
end
println("\n")
group_number+=1
end

Group 1 :
Alice, Bob, Charles, Evan,

Group 2 :
Charles, Dick, Ron, Alice,

Group 3 :
Evan, Dick, Ingrid,

Let's re-run this for the example we constructed by hand above, and compare it against our earlier model

Just to recap, here's the earlier example :

Group 1 (Alice, Bob and Charles)

Alice owes Bob \$10 Bob owes Charles \$10

Group 2 (Dick and Charles)

Charles owes Dick \$10 Let's set up the$T$and disallowed_transanction matrices manually. In [31]: T = [ 0 10 0 0; 0 0 10 0; 0 0 0 10; 0 0 0 0; ]; unique_users = ["Alice", "Bob", "Charles", "Dick"] n = 4 # Dick can't transanct with either Bob or Alice disallowed_transanctions = [ (1,4) (4,1) (2,4) (4,2) ];  In [32]: S = find_S_matrix_grouped(T, n, unique_users) for i in unique_users for j in unique_users if S[i,j] != 0 println("$i owes $j \$$(S[i,j])") end end end  Optimal Alice owes Charles 10.0 Charles owes Dick 10.0  This is what we expected, but now let's just compare it to what we would get if we just used our earlier solver. In [33]: S = find_S_matrix(T, n, unique_users) for i in unique_users for j in unique_users if S[i,j] != 0 println("i owes j \$$(S[i,j])") end end end  Optimal Alice owes Dick$10.0

Of course, it doesn't know anything about groups, so it just (wrongly) suggests Alice should pay Dick.

### 3.3 Determining the Upper Bound¶

Till now, we've been using an arbitrarily large value for the upper bound on the problem, but we can actually determine this by just looking at the raw transanctions.

The upper bound in our case is essentially the largest single transanction that could occur in the settlement phase.

This can also be thought of as the maximum amount of money that a single person has to receive. We don't have to worry about the maximum amount of money a single person owes, because even if it's greater than the maximum someone has to receive, that particular value will never occur in the settlement phase because the amount owed would be split between at least two people.

In [34]:
function find_upper_bound(T, n)
colsums = []
for i in 1:n
colsum = sum(T[j,i] for j in 1:n)
push!(colsums, colsum)
end
return maximum(colsums)
end

Out[34]:
find_upper_bound (generic function with 1 method)
In [35]:
T, n, unique_users = get_T_from_csv("dense.csv");

In [36]:
find_upper_bound(T, n)

Out[36]:
1190.0

It's trivial to integrate this reasonining into our function that actually solves the problem, so for brevity, it's not done here.

## 4. Results and Future Work¶

So we looked at how we can use integer programming in order to simplify debts across multiple groups. This can make a big difference in the number of transnactions required to get to a settled state.

There's many different aspects of this problem that can be explored further :

• Are there other, simpler methods of modelling this problem? The number of variables we use is of the order $n^2$ for every $n$ users. Therefore this approach will not scale to hundreds or thousands of transanctions. I had initially explored the possibility of modelling this as a flow problem, which would simplify the model greatly. I think this is something definitely worth purusing.
• Zero diagonal. One of the characteristics of the problem we could make use of is the fact that the main diagonal of both the $T$ and $S$ matrix have to always be zero. This would require changing the model quite a bit, but it would reduce the problem to requiring n(n-1) variables, which on paper doesn't seem like a big difference, but it can make a difference for large values for $n$.

## 5. Interactive Mode¶

Do you have debts you need simplified? Go ahead and put your transanctions in and run the cell after that to find the quickest way to get to a settled state with your friends!

In [37]:
my_transanctions = [
# format : "owed by,owed to,amount",
# don't forget the , at the end
# sample :
"Tony,Bob,10",
];

In [38]:
function get_T_manual()
transanctions = []
all_users = []
lines = my_transanctions
for line in lines
owed_by, owed_to, amount = split(line, ",")
push!(transanctions, [owed_by, owed_to, parse(Int, amount)])
push!(all_users, owed_by)
push!(all_users, owed_to)
end
unique_users = sort(unique(all_users))
n = length(unique_users)
T = NamedArray(zeros(n,n), (unique_users, unique_users), ("OWED BY", "OWED TO"));
for transanction in transanctions
owed_by, owed_to, amount = transanction
T[owed_by, owed_to] = amount
end
return T, n, unique_users
end
T, n, unique_users = get_T_manual()
S = find_S_matrix(T, n, unique_users)
for i in unique_users
for j in unique_users
if S[i,j] != 0
println("$i owes$j \(S[i,j])")
#println("\"$i,$j,$(S[i,j])\",") end end end  Optimal Tony owes Bob$10.0

format: plain


text:
• # question_answer 75)   Explain the non-linear shape of ${{H}_{2}}S$ and non-planar shape of $PC{{l}_{3}}$ using valence shell electron pair repulsion theory.

Answer:

The Lewis structure of ${{H}_{2}}S$ is : S-atom is surrounded by four electron pairs (two bonded and two lone pairs). These four electron pairs adopt tetrahedral arrangement. The presence of two lone pairs brings distortion in the molecule on account of repulsion with bonded pairs of electrons. Thus, the shape of ${{H}_{2}}S$ molecule is V-shaped and not linear. The Lewis structure of $PC{{l}_{3}}$ is : P-atom is surrounded by four electron pairs (3 bonded and one lone pair). These four pairs adopt a tetrahedral geometry. Due to the presence of lone pair, $PC{{l}_{3}}$ has a distorted tetrahedral geometry. Thus, it is pyramidal in shape and not non-planar shape.

You need to login to perform this action.
You will be redirected in 3 sec


text:
Properties of quadratic equations and inequalities

Question Sample Titled 'Properties of quadratic equations and inequalities'

If the quadratic equation ${h}{\left({x}+{5}\right)}^{{2}}={k}$ , where ${h}$ and ${k}$ are constants, has distinct real roots, then which of the following may be true?

 I. ${h}$ $>{0}$ and ${k}>{0}$ II. ${h}$ $<{0}$ and ${k}<{0}$ III. ${h}$ $>{0}$ and ${k}={0}$

A
I and II only
B
I only
C
I, II and III
D
II and III only

 ${h}{\left({x}+{5}\right)}^{{2}}$ $={k}$ ${h}{\left({x}^{{2}}+{10}{x}+{25}\right)}$ $={k}$ ${h}{x}^{{2}}+{10}{h}{x}+{25}{h}-{k}$ $={0}$ ∵   The quadratic equation has distinct real roots. ∴  $\Delta>{0}$ ${\left({10}{h}\right)}^{{2}}-{4}{\left({h}\right)}{\left({25}{h}-{k}\right)}$ $>{0}$ ${100}{h}^{{2}}-{100}{h}^{{2}}+{4}{h}{k}$ $>{0}$ ${h}{k}$ $>{0}$ ∵  ${h}{k}>{0}$ ∴  ${\left({h}>{0}\right.}$ and ${k}>{0}{)}$ or ${\left({h}<{0}\right.}$ and ${k}<{0}{)}$ , and also ${k}\ne{0}$ .

專業備試計劃

Level 4+ 保證及 5** 獎賞

ePractice 會以電郵、Whatsapp 及電話提醒練習

ePractice 會定期提供溫習建議

Level 5** 獎勵：會員如在 DSE 取得數學 Level 5** ，將獲贈一套飛往英國、美國或者加拿大的來回機票，唯會員須在最少 180 日內每天在平台上答對 3 題 MCQ。

Level 4 以下賠償：會員如在 DSE 未能達到數學 Level 4 ，我們將會全額退回所有會費，唯會員須在最少 180 日內每天在平台上答對 3 題 MCQ。

FAQ

ePractice 是甚麼？

ePractice 是一個專為中四至中六而設的網站應用程式，旨為協助學生高效地預備 DSE 數學（必修部分）考試。由於 ePractice 是網站應用程式，因此無論使用任何裝置、平台，都可以在瀏覽器開啟使用。更多詳情請到簡介頁面。

ePractice 可以取代傳統補習嗎？

1. 會員服務期少於兩個月；或
2. 交易額少於 HK\$100。

Initiating...

HKDSE 數學試題練習平台


text:
# A barometer kept in an elevator reads 76cm when it is at rest. If the elevator goes up with

Q: A barometer kept in an elevator reads 76 cm when it is at rest. If the elevator goes up with increasing speed, the reading will be
(a) zero

(b) 76 cm

(c) < 76 cm

(d) > 76 cm

Ans: (c)

Sol: As atmospheric Pressure is same in both the condition

$\displaystyle \rho g h = \rho (g + a)h’$

$\displaystyle h ‘ = \frac{g }{g + a} h$

Hence , h’ < 76 cm


text:
## Simultaneous stabilization of linear systems

### Vincent D. Blondel

Springer Verlag, Heidelberg, 1994
Lecture Notes in Control and Information Sciences, ISBN: 3-540-19862-8

Review in Zentralblatt MATH (795.93083):

"The central question investigated in this beautiful research monograph is the problem of simultaneous stabilization: Let p_1, p_2, ..., p_k be k linear systems. Under what conditions is it possible to find a controller C which is stabilizing for each p_i (i= 1,... , k).
The traditional examples of the k systems include nominal systems, with its failed modes, an airplane which has natural operating points (it takes off, flies and lands) etc. Unexplored but potentially fruitful are systems of dynamics of different national incomes, p_1, p_2, ..., p_k. Can the U.N. -- World Bank provide a stabilizing controller C for each p_i.

The author provides the current state of attempts to solve this problem. Simultaneous stabilizability question for k=2 is reported to be fully solved. There exists no tractable condition for simultaneous stabilization of three or more linear systems. Thus if we allow only finite combination of rational, logical and sign test operations on the coefficient of the k (k > 2)$systems, it is not possible to find necessary and sufficient conditions for simultaneous stabilization. The text is very lucid and carefully written with an excellent summary and bibliography at the end of each chapters. It contains a mathematical appendix. The English'' prologue and Epilogue is beautiful. It succeeded in its attempt to explain mathematical theory to the first man in the street. This research monograph is recommended to control theorists, applied mathematicians and systems engineers." From the introduction: "This monograph focuses on the question of simultaneous stabilization of scalar linear systems: Let$p\sb 1,p\sb 2,\cdots, p\sb k$be$k$scalar linear systems. Under what condition does there exist a controller$c$that is stabilizing for each$p\sb i (i=1,\cdots,k)$? We are most interested in the existence question and we make only a few remarks on constructive procedures. The volume consists of 6 chapters and 3 appendices. The chapters are as independent of each other as possible. They all start with a short introduction and end with a concluding section that summarizes the content and gives bibliographical references. Chapter 2: Stabilization. This chapter contains all the basic definitions: stabilization, strong, bistable and simultaneous stabilization, etc. There exist several conceptual ways to look at stabilization and we present two of these: the algebraic factorization approach and the geometrical avoidance interpretation. Chapter 3: Youla-Kucera parametrization. The Youla-Kucera parametrization of the set of all stabilizing controllers of a given system is the stepping stone to many control design strategies. It is presented here and used to derive equivalences between strong, bistable and simultaneous stabilization. Chapter 4: Necessary conditions. A rational function that has no poles in the right half plane has no poles on the positive real axis. This elementary observation is the guideline used in the chapter to obtain necessary conditions for simultaneous stabilization. If all the closed-loop transfer functions associated to the systems$p\sb i (i=1,\cdots,k)$in feedback with a controller$c$have no poles in the right half plane, then this controller is such that all the closed loop transfer functions have no poles on the positive real axis. Such a controller is said to be$\bold R\sb {+\infty}$-stabilizing for$p\sb i (i=1,\cdots,k)$. The existence of an$\bold R\sb {+\infty}$-stabilizing controller is thus a necessary condition for the existence of a stabilizing controller. In this chapter we derive tractable necessary and sufficient conditions for simultaneous$\bold R\sb {+\infty}$-stabilizability of$k$systems and for strong and bistable$\bold R\sb {+\infty}$-stabilization. These conditions encompass all known necessary conditions and are expressed in the form of interlacement properties on poles and zeros of the systems and of related rational functions. Chapter 5: Sufficient conditions. In this chapter we show that the necessary conditions derived in Chapter 4 are also sufficient for simultaneous stabilization of two systems and for strong stabilization. Then we show that this property does not flow on to bistable stabilization or to simultaneous stabilization of more than two systems: three systems that are simultaneously$\bold R\sb {+\infty}$-stabilizable are not guaranteed simultaneous stabilizable. The conditions given in Chapter 4 for bistable or simultaneous stabilization of$3$or more systems are not sufficient. The proof of this result makes crucial use of results from analytic function theory. In the last two sections we provide examples of sufficient conditions that are related to the avoidance concept given in Chapter 2 and to$H\sb \infty\$ control design.

Chapter 6: Necessary and sufficient conditions. The result contained in this closing chapter is that the simultaneous stabilization question is rationally undecidable: it is not possible to find necessary and sufficient conditions for simultaneous stabilization of three or more systems that involve only a finite combination of rational operations (additions, subtractions, multiplications and divisions), logical operations (and' and or') and sign test operations (equal to, greater than, greater than or equal to, etc.) on the coefficients of the three systems."


text:
# orthogonal vectors point of intersection

• Dec 23rd 2009, 10:24 AM
thesys
orthogonal vectors point of intersection
Hello,

Hello,
I have two points p and q. i need to find point r which is defined by the intersection of the vectors pr and rq. the vectors are orthogonal, so that <pr,qr> = 0

how can i get the point r ? How can I define my lines so i can make pl1 = ql2
• Dec 23rd 2009, 11:32 AM
Soroban
Hello, thesys!

If I understand the problem, there is an infinite number of answers.

Quote:

I have two points $P$ and $Q.$
i need to find point $R$, defined by the intersection of $\overrightarrow{PR}$ and $\overrightarrow{QR}.$

The vectors are orthogonal, so that: . $\overrightarrow{PR}\cdot\overrightarrow{QR} \:=\:0$

How can i get the point $R$ ?

Code:

              * * *    R           *          o         *          *    *       *      *        *             *       *  *                *     P o - - - - * - - - - o Q

We have a semicircle with diameter $PQ.$

Let $R$ be any point on the semicircle.
Draw chords $PR$ and $QR.$

Then: . $PR \perp QR$

• Dec 23rd 2009, 11:55 AM
thesys
Quote:

Originally Posted by Soroban
Hello, thesys!

If I understand the problem, there is an infinite number of answers.

Code:

              * * *    R           *          o         *          *    *       *      *        *             *       *  *                *     P o - - - - * - - - - o Q

We have a semicircle with diameter $PQ.$

Let $R$ be any point on the semicircle.
Draw chords $PR$ and $QR.$

Then: . $PR \perp QR$


text:
How to use memory without a trace

Harry Buhrman is a quantum complexity theorist. That is, if you take a measure of him, you will sometimes find that he is doing complexity theory, sometimes quantum computing, sometimes both, and sometimes neither. Whichever, he has been doing it great for a quarter century. He has a joint paper in STOC 2014, which begins today, that is part complexity and part—something else.

Today Ken and I want to talk about the beautiful new concept and results in that paper. They can be viewed as complexity theory results: this model of computation can compute such-and-such. They can also be viewed as security results: malware programs can do things to avoid detection.

We are not suggesting that Harry and his co-authors—Richard Cleve, Michal Koucký, Bruno Loff, and Florian Speelman—are up to something nefarious. Their paper is called Computing with a full memory: Catalytic space. No hint of security or anything diabolical here. Indeed they seem to think of the result more as a result about chemical reactions, hence the term “catalytic.”

We beg to see their result in a different light. Of course, the authors can call it whatever they like: one of the perks of writing the paper. Of course, we can interpret the result also in any way that we like: fair is fair. Seeing it as related to malware and security is meant as a compliment, and we certainly hope that Harry and his co-authors take it that way. Or perhaps they might take it in a complementary way.

## A New Model of Computation

So what is their new concept? It is a new type of model of computation. We quote them:

Imagine the following scenario. You want to perform a computation that requires more memory than you currently have available on your computer. One way of dealing with this problem is by installing a new hard drive. As it turns out you have a hard drive but it is full with data, pictures, movies, files, etc. You don’t need to access that data at the moment but you also don’t want to erase it. Can you use the hard drive for your computation, possibly altering its contents temporarily, guaranteeing that when the computation is completed, the hard drive is back in its original state with all the data intact?

They go on to qualify an important point:

One natural approach is to compress the data on the hard disk as much as possible, use the freed-up space for your computation and finally uncompress the data, restoring it to its original setting. But suppose that the data is not compressible. In other words, your scheme has to always work no matter the contents of the hard drive. Can you still make good use of this additional space?

The short answer is yes. In order to state their result we need to formalize their model a bit. Their model leans on a restricted class of ordinary random-access machine (RAM) programs with the standard numerical ${+,-,*}$ operations, and registers ${r_1,r_2,r_3,\dots}$ Their transparent programs have instructions only of the form

$\displaystyle r_i \text{ \tt += } r_j*r_k \qquad\text{and its inverse}\qquad r_i \text{ \tt -= } r_j*r_k,$

where ${j,k \neq i}$. If we restrict ${i > n}$ for some ${n}$, then the initial contents ${x_1,\dots,x_n}$ of ${r_1,\dots,r_n}$ are read-only inputs. The main point is that the model also allows ${r_{n+1}}$ through the last register ${r_N}$ to have arbitrary initial values. A program ${P}$ computes ${y = f(x_1,\dots,x_n)}$ if the initial value ${r_N}$ is transformed to ${r_N + y}$.

The stealthy aspect is that if you allocate one more register ${r_{N+1}}$, then you can write the program

$\displaystyle P' = r_{N+1} \text{ \tt -= } r_N;\;\;P;\;\;r_{N+1} \text{ \tt += } r_N;\;\;P^{-1},$

where ${P^{-1}}$ applies the inverse of each instruction in ${P}$ in reverse order. This computes ${y}$ into ${r_{N+1}}$ while leaving all previous registers unaltered. The trick works even if the values belong to a finite ring. With a finite ring we can model the registers by regularly-spaced finite segments of a tape, and thus hope to control the programs by Turing machines.

## The Idea

With ${s(n)}$ tape cells, one can address any location on a tape of size ${2^{s(n)}}$. Thus a log-space Turing machine can control a polynomial number of registers, and also read the input tape. If the programs ${P_n}$ to run on inputs of length ${n}$ are uniform enough to generate each successive instruction in logspace, then our Turing machine can execute them. Our machine can also save an ${O(\log n)}$-sized portion of a larger tape, introduce its own value there, and restore the original later.

Now suppose the input holds the adjacency matrix ${A}$ of a directed graph on ${n = 2^k}$ nodes. If we can square ${A}$ iteratively ${k}$ times, then we can tell whether there is a path from node ${1}$ to node ${n}$, which solves a problem complete for nondeterministic logspace. Squaring ${A}$ has very regular polynomial-length transparent programs.

The final hitch is that the final entry ${A^n(1,n)}$, and entries along the way, can have values up to ${n^n}$. One solution is to compute modulo ${p}$ for a long enough logspace-uniform sequence of small primes ${p}$; then by the Chinese remainder theorem we can accept iff any one of them gives a non-zero value mod ${p}$. Since ${p}$ has ${O(\log n)}$ length our machine can preserve the original (arbitrary) contents ${r_N}$ of the destination register, recognize that the value was changed, and restore ${r_N}$.

## Main Results

Their formal model thus augments the standard Turing machine model—which has work tapes and read-only input—with an auxiliary tape. They call the machine catalytic if it computes the same value for every possible initial setting of the auxiliary tape, and that tape is unchanged at the end. They define ${\mathsf{CSPACE}(s(n))}$ to be the class of sets computed by catalytic Turing machines whose work-tape is bounded by ${s(n)}$ tape cells, and whose auxiliary space is bounded by ${2^{O(s(n))}}$ cells. Write ${\mathsf{CL}}$ for ${\mathsf{CSPACE}(\log n)}$. This is a deterministic log-space class, yet by the above argument, it contains non-deterministic logspace: ${\mathsf{NL \subseteq CL}}$. In fact, they prove something stronger:

Theorem 1 ${\mathsf{TC^1 \subseteq CL \subseteq ZPP}}$.

Here ${\mathsf{TC^1}}$ is the class of languages accepted by uniform log-depth threshold circuits, and ${\mathsf{ZPP}}$ is the class of languages with programs that halt in expected polynomial time and are always correct when they halt. The former inclusion builds on wonderful known connections between arithmetical and Boolean complexity classes using counting as a go-between, and their paper gives a great survey of them. The latter inclusion exploits the fact that since the machine is deterministic, it can’t have any overlap in configuration between computations with different initial settings of the auxiliary tape, because the flowing together would make the final auxiliary tape contents wrong for one of them. Since the remaining space is logarithmic, the average number of disposable configurations per setting is polynomial. Hence a randomized TM that sets the auxiliary tape randomly runs in expected polynomial time. They also prove:

Theorem 2

• If the Exponential Time Hypothesis is true, that is ${\mathsf{SAT} \notin \mathsf{BPTIME(2^{o(n)})}}$, then ${\mathsf{SAT} \notin \mathsf{CSPACE(o(n))}}$.

• There is an oracle ${A}$ such that ${\mathsf{CL}^A = \mathsf{PSPACE}^A}$, indeed ${\mathsf{CSPACE}^A(s(n)) = \mathsf{DSPACE}^A(2^{O(s(n))})}$ for every ${s(n) \geq \log n}$.

• There is an oracle ${B}$ such that ${\mathsf{NL}^B \not\subseteq \mathsf{CL}^B}$.

The last statement seems jarring since above we proved ${\mathsf{NL \subseteq CL}}$, but it plays on known issues with defining oracles for space-bounded machines. That’s the very last result in their paper, but…

## More Applications?

Note what they show is that memory that is already in use can be used for temporary storage during a computation. This is why we see their result as related to security and malware.

Imagine the following scenario. You are a malware program. You do not want to ask for additional storage, since this could alert the system that you are running. Instead you plan to use memory that is currently being used for storage of objects like movies or whatever. The key is that you use this memory to expand your power of computation, but when done you leave the memory in the exact state that it was in when you started.

Clearly, this a great for a malware program that needs lots of memory during its execution. It can avoid making any system call to allocate additional memory, yet can use tons of memory while it is executing. This stealth type of behavior seems quite cool, quite powerful, surprising, and if practical could be used in real attacks on systems.

## Open Problems

As they say, how does ${\mathsf{P}}$ relate to ${\mathsf{CL}}$? As we say, if every polynomial-time program could be “catalyzed” by a program with a small own-footprint, what greater dangers would we face?

1. June 2, 2014 6:35 am

Nothing diabolical? I can think of one demon who might be interested – Maxwell’s. A bank of memory full of zeroes is a useful resource for such a demon. This seems to suggest that memory full of random stuff (equivalent to a system at thermodynamic equilibrium) is just as good. Does this mean that we can avoid the heat death of the universe? Where’s the catch?

• June 2, 2014 11:03 am

That’s a neat idea, and spurs me to read Roger Penrose’s Cycles of Time which I hear talks about creatio ex equilibrio.

June 2, 2014 7:23 am

In future versions of Linux, the home and the swap partitions will be one and the same…