Dynamic Model for Population Distribution and Optimum Immigration and Job Creation Policies

In this paper we demonstrate that by use of modern Systems and Optimal Control theory, it is possible to formulate optimum immigration and job creation strategies while maintaining population level close to certain pre-specified targets. With this objective in mind, we consider a simplified dynamic model based on a previous model developed in (Ahmed and Rahim, 2001:325-358) to describe the population distribution in Canada. Numerical results demonstrate that the model population is in close agreement with the actual population. This model was then used to formulate a control problem with immigration and job creation rates being the decision (control) variables. Using optimal control theory, optimum immigration and job creation policies were determined. Results are illustrated by numerical simulation and they are found to be very encouraging.


Introduction
The main objective of this paper is to demonstrate that by use of modern Systems and Optimal Control theory, it is possible to formulate optimum immigration and job creation strategies while maintaining population level close to certain pre-specified levels.It is reasonable to think that the immigration should be tied with the demand for manpower and the availability of jobs and ability to create new jobs.
Since the very first demographic model was proposed by Malthus (1798), there has been a great deal of research on building population models (Sharpe and Lotka 1911: 435-438;Pollard, 1973:23-26;Das Gupta, 1978: 367-379;and Lee, 1974: 563-585).Most of this demographic research has concentrated on fertility and birth function and from it predict population change.This is fundamental and it is required in modeling any population dynamics.In contrast, our research will focus on constructing dynamic models for population distribution in Canada by following the systems and control theoretic approach as developed in (Ahmed and Rahim, 2001: 325-358).The methodology proposed in this paper may be useful for the Government of Canada to formulate optimum immigration and job creation policies.

Dynamic Model for Population Distribution and Optimum Immigration and Job Creation Policies
In general, population changes with time because of birth, death, immigration and emigration etc.If the total population is divided into several age groups, there is a continuous process of transition from one age group to the next.Following Statistics Canada, total population is divided into three age groups.Group one, denoted by G1, includes population in the age group (0,14], group G2 includes population in the age group [15,64], and finally group G3 absorbs the reaming population.To capture the temporal variation of population in each group one must build a dynamic model.Such a model can be constructed if basic parameters like birth rate, death rates, and transition rates from one age group to the next are available.Very often these parameters are not readily available; only the population data in each individual age group is available for a given period of time.Based on this information, estimation and identification of unknown parameters are explored using the mathematical model mentioned above (Ahmed and Rahim, 2001: 325-358).
Out of all the parameters mentioned, only immigration rate can be easily monitored and controlled by the Canadian government (Beaujot and Matthews 2000).(By immigration rate of any age group G1, G2, G3, we mean the intake of new immigrants of that age group per week per unit of population of the same group).Immigration plays an important role in avoiding population decline and maintaining the necessary labor force for economic growth (Beaujot, 2002).According to current practice, each year the target level of immigration, for the following year, is proposed.Obvious question is, how is the target level of immigration determined?In general, it is influenced by changes in manpower requirement and other socio-economic and political conditions of the country.If there are not enough job positions to match the level of immigration, it may cause many social and political problems.Currently, immigration and unemployment in Canada is an interesting and important topic (Siklos and Marr, 1998:127-147).Some papers have been published investigating the immigration and UI system (Hornstein and Yuan 1998) as well as immigration and the rate of growth of population and labor force (Denton, Feaver and Spencer 1997).In our research, we consider immigration rate and job creation rate as the control or decision variables which the federal government, in consultation with its provincial counterparts, can adjust as required to change the dynamics of population distribution.Here, by job creation rate we mean the number of new jobs created during a week divided by unemployed population in age group G2 during the same week.
The objectives of this paper are: (1) to construct a dynamic model to describe the population distribution of Canada in the three age groups as mentioned above, (2) to identify the unknown parameters by minimizing the identification error (defined later) when the parameters are not available, (3) to optimize the immigration policy (rate) by minimizing an objective functional called the (cost function), and (4) to optimize immigration and job creation rates by minimizing a similar cost function (described later).
The paper is organized as follows.Section 2 presents a population model.Section 3 applies the model to identify the unknown parameters by minimizing the error between the model response and actual population.Section 4 explores the optimum immigration policy using the population model (constructed).Section 5 includes in the population model the unemployment factor and seeks for optimum immigration and job creation rates.

Dynamic Model
We use similar models as developed by (Ahmed and Rahim, 2001: 325-358).In that paper, they presented a set of (dynamic) mathematical models to represent the population distribution in general.The dynamic models presented in this paper are modified according to the actual population grouping and characterization used by Statistics Canada.This is then used to describe the population variation (dynamics) in Canada.A question that was raised by the reviewers of this paper is why only 3 age groups?In our original model (Ahmed and Rahim, 2001: 325-358), as mentioned above, we proposed 4 age groups with G1 (0, 13) consisting of children below age 13, G2 [13,19) consisting of teenagers between 13 to 19 and G3 [19, 65) comprised of adults from 19 to 65. Population above 65 is grouped as G4.This grouping allows the federal as well as the provincial governments to estimate the cost of administering various government services directed to various age groups and thereby formulate reliable policies and allocate resources as required.In order to be able to compare our results with the available data from Stat-Canada, we regrouped the population according to Stat-Canada practice reducing from 4 age groups to 3. In fact our modeling procedure allows segmentation of the population into as many groups as required without any difficulty.

Population Dynamics
Following Statistics Canada, the Canadian population is divided into three age groups.Group G1 consists of all children up to age 14, and we let the population in that group at any time t be represented by x .Group G2 consists of the population between ages 15 to 64, and we denote the population in that age group at any time t by ) ( 2 t x .Population over age 65 is grouped as G3, and the population count in this age group is represented by The reason population is partitioned into these different age groups is that many Government programs (for example; child care, education, health, old age pension, unemployment insurance, welfare etc.) and the cost of administering those programs and services are strongly dependent on the population distribution.

Mathematical Model
The growth rate of population of group G1 is given by where b denotes the birth rate due to population of age group G2.We assume that contribution to birth rate due to population below age 15 and above age 65 is negligible.The parameter 12 denotes the maturity rate or transition rate from age group G1 to age group G2 (see section 3 for exact definition), d1 denotes the child mortality rate, e1 and i1 denote the child emigration and immigration rates respectively.Similarly the growth rate of G2 can be written as where the parameter 23 denotes the transition rate from age group G2 to age group G3, d2 denotes the mortality rate of the age group G2, e2 and i2 the corresponding emigration and immigration rates respectively.The growth rate of G3 can be written as where the parameters d 3 , e3 , i3 are, respectively, the mortality rate, emigration rate, and immigration rate of population in age group G3.We can write the above population model in the Canonical form as follows: . , 0 )) ( ( (5) denotes the population vector in the three age groups.The vector function denotes the vector field representing the population growth rates.

System Based on Available Parameters
Yearly data on birth, death, immigration, emigration and population of all the three age groups for the period 1972 to 2002 was obtained from Statistics Canada.This data was then converted to weekly values.Here, 12 stands for the weekly transition rate from G1 to G2 and 23 from G2 to G3.The coefficient b stands for the weekly birth rate, stand for weekly mortality rates, stand for weekly immigration rates, and e e e 3 2 1 , , weekly emigration rates for the groups G1, G2, and G3, respectively.For example consider the rate 12 for the week wi , i = 1, 2,… , 52 for any given year.This is given by ) ( 1 where G1(14, 1yr) denotes the number of 14 year old members of the population G1 during the year and G1(i) denotes the population of group G1 during the ith week of the year.Consider another rate like the death rate d2 .This is given by ) ( 2 where ) 1 ( 2 yr d denotes the number of deaths in age group G2 during the year under consideration and G2(i) stands for population of group G2 during the i -th week of the same year.All the remaining rates are determined in the same way.

Numerical Results Based on Available Parameters
Comparison of the actual population with the model population (over the years , using the infinitesimal rates computed from the data obtained from Statistic Canada, is shown in Figures 1-3.The dashed lines represent the actual population and solid lines represent the model population.
Based on the assumptions mentioned above, we conclude from the numerical results that the model population is fairly close to the actual population though it falls short of exact match.The mismatch is possibly due to the presence of no linearity of true population model and numerical errors in computing the rates from the raw data.Using neural networks it is possible to construct more accurate nonlinear models.We leave it for future investigation.

Parameter Identification
In the previous section, parameters, like birth rates, death rates and transition rates, were assumed to be given or computed approximately from available population data.But in some cases, the actual parameters mentioned above may not be available except the population data for a particular period of time.In that situation, the methodology of system identification provides an effective tool for estimating these unknown parameters (Ahmed, 1976).They can be determined also by an alternative approach provided by optimal control (or decision) theory (Ahmed, 1988) and this is what we use in this paper.The term "control" is generally used in engineering literature and currently it is also widely used in mathematical sciences.However the concept is much broader than what is implied by physical sciences.It includes actions of any kind (for example, passing of new legislations, introduction of city bylaws, changes in prime lending rates, bank interest rates, diversion of high ways and rivers, deforestation, plantation etc.) that may impact the society in anyform.Optimal control theory provides a rigorous mathematical tool for determination of policies that produce the best result within the limits of social and financial resources.

Problem Formulation and its Solution
We formulate the identification problem as an optimal control problem.The unknown parameters are found by minimizing the error between the model response and the actual data.Population in the three age groups are known for a period of time, say I [0, T], including the weekly emigration and immigration rates.We need to estimate the weekly birth and death rates for the three age groups and the transition rates from G1 to G2 and from G2 to G3.These six parameters, which are functions of time, are to be identified (estimated).The state equation, as described before, is given by with the unknown parameters denoted by The objective is to find the vector that minimizes the error ) ( J .We use optimal control theory (Pontryagin Minimum Principle) (Ahmed, 1988: 258- where the adjoint state (co-state) (also known as Lagrange multiplier) is given by the solution of the adjoint system, written explicitly as, subject to the terminal condition given by, (t)= 0, that is, The gradient vector is given by Using the state and co-state equations and the gradient as defined above, one can write an algorithm for numerical solution of the optimization problem.There are three different algorithms as suggested in (Ahmed, 1988: 302-304) any one of which can be used to compute the best parameter 0 (minimizing the functional ( 9)).For a complete description of these algorithms see [ (Ahmed, 1988: 302-304) It is known from basic calculus that for minimizing or maximizing a functional of certain decision variables, one equates the gradient (first derivative) to zero and solves the resulting (generally) nonlinear equation for the unknown variable.If the unknown variable is a vector z of finite dimension this is pretty simple.However if there are also some side constraints that the unknown variable is required to satisfy, even the finite dimensional problem may present sufficient difficulties.In this case one introduces the so called Lagrange multiplier to convert the constraint problem into an unconstraint problem which is then solved for the unknown z and the Lagrange multiplier simultaneously to find the optimal z0 .In other words the Lagrange multiplier helps in determining the best direction of approach towards the optimum.Optimal control theory is a far reaching generalization of this basic philosophy which is applicable to problems involving infinitely many variables (or continuum of variables) subject to dynamic and static constraints.Loosely speaking, optimal control theory presents a systematic and powerful tool for computing the gradient successively thereby providing the direction of steepest decent for minimization (steepest accent for maximization) of any objective functional (subject to dynamic and static constraints) leading to optimum controls or decisions.All the optimization problems considered in this paper involve decision variables which are functions of time and hence infinitely many.Optimal control theory provides the necessary gradient through the Hamiltonian H and the co-state variable (also called the Lagrange multiplier).This is the basic tool used here throughout the paper.For more details the reader may see [Ahmed, 1988:chapter 6, p232].

Numerical Results Based on Identified Parameters
For the three age groups, Figures 4-6 show the actual and model population, the later based on identified parameters (birth rate, death rates and transition rates).
Because of the scale it appears that there is no difference between the actual and the model population.To see the actual difference in these figures, a segment of Fig 4 covering week 700-720 is presented in Table 2.In fact the maximum error over the entire period of (31years) is about ± 327, which is considered to be an excellent match.

Comparison of Results Based on Available and Identified Parameters
The numerical results based on actual parameters are shown in Figures 1-3, and those based on identified parameters are displayed in Figures 4-6.It is clear that the model population based on identified parameters matches the real population much better.

Optimum Immigration Policy
Some developing countries like China and India are over populated.This causes many social problems in areas like education, employment, transportation and environment.In order to pursue better living standards, people from these developing countries are eager to migrate to developed countries.On the other hand, some developed countries, like Canada, Australia and New Zealand accept immigrants from other countries to prevent their own population decline and maintain a labor force for economic growth (Beaujot, 2002).Canada has some very clear immigration policies.In the fall of each year, the immigration level for the coming year is set by the Canadian Government in consultation with their provincial counterparts.Using systems and control theory as proposed here one can determine the optimum immigration levels subject to any number of constraints that may be applicable.This is formulated in the following subsection.

Problem Formulation and its Solution
Suppose our objective is to maintain the total population above a minimum level and not exceed far above a maximum level.Ideally the immigration policy should be one that maintains the population variation within the specified band determined by the minimum and the maximum levels.This is to be achieved by adopting an appropriate immigration rate (policy) for the age group G2 by minimizing a cost function, which penalizes whenever the population level is outside the boundary of the target set as described above.Whenever the total population is within the lower and upper limits of the target set [ ], the cost function given by the expression ( 18) is zero.When it is below (above) the lower boundary (the upper boundary), a quadratic penalty is imposed on the mismatch with weight ) ( 1 2 as seen in the expression (18).These parameters can be chosen by the planner assigning priorities as required.In case the population is below the lower boundary, to avoid population decline one may choose larger values for 2 relative to 1 .On the other hand if the population is above the upper boundary of the target set, to avoid explosive growth and associated social pressure one may choose larger values for 1 relative to 2 .
Parameters, such as weekly birth and death rates; transition rates from group G1 to group G2; and group G2 to group G3; emigration rates; and the lower and upper boundaries of the target population, are given.The immigration rate u of the age group G2 (adult population) is the control variable to be determined.
The state equation, including immigration rate as the control variable, is given by where u denotes the (weekly) immigration rate.The parameters p 1 and p 2 are nonnegative fractions determining the number of accompanying children and seniors.The total number of new immigrants per unit of time is given by the sum of adults and their accompanying children and seniors, which equals In view of the objective as described above, the cost function is defined as follows: The corresponding adjoint system is given by or equivalently, where, for convenience of presentation, we have introduced the function as given by ).
Since the cost function does not contain any terminal constraint, the terminal condition for the co-state ( 23) is given by (T) = 0.In this case the gradient vector is given by

Dynamic Model for Population Distribution and Optimum Immigration and Job Creation Policies
Using the state and co-state equations along with the gradient as defined above, we use the same algorithm as indicated in section 4.1, to determine optimum policies.

Numerical Analysis and Results
Without Specified Target (Fig. 7-Fig.9) In this subsection we present numerical results on population growth corresponding to actual immigration and emigration data (Fig. 7) obtained from Stat.Canada.Here no target population was specified.Figure 8 shows the actual population dynamics in Canada from 1976 to 2002 as reported by Stat.Canada.Figure 9, describes the corresponding model population based on identified birth, death, and transition rates, including given immigration and emigration rates as shown in Fig. 7A-7B.The results are quite close.Clearly Fig. 9 shows the natural population dynamics without any external intervention (that is without any other control).Thus according to current immigration policy, the Canadian population would monotonically increase with time.
With Specified Target (Fig. 10-Fig.14) In this subsection we present numerical results corresponding to population targets as described in section 5.1.Using the system (16) with immigration rate u as the control variable and (18) as the objective functional, and using the optimization procedure as described in section 5.1, we obtain the optimum immigration policy.The results corresponding to the optimum immigration policy are displayed in Figures 10-14.
(A) Fixed Target  Optimum control (or immigration) policies corresponding to the target sets ] are shown in figures 10 and 11.The curves 10B and 11B are the expanded versions of 10A and 11A respectively showing the detailed transition from high to low immigration rates.Note that the immigration rate corresponding to the target set T 2 remains at its maximum admissible value for a longer period of time compared to that of target set T1 because the lower and the upper boundaries of the set T2 are above those of T1 .It appears from this result that optimum immigration policy is to keep the immigration rate (actual number = rate x2 ) at its highest admissible level during the early period of the planning horizon and then rapidly reduce to zero.This translates into maximum 500,000 (approximately) annually during the early period.In control theory this kind of phenomenon is known as bang-bang control.This shows that once the lower limit of the target set is met there is no reason for further immigration unless humanitarian concerns are added to our cost function.This is entirely due to the population goal and the maximum immigration rate used in our example.If there is no specified target, we have already seen at the beginning of this subsection that population will continue to increase with the increase of immigration rate.
The graphs of Figure 12 show the corresponding population growth for target sets T1 and T2 .The cost function ( 18) is designed so as to reach the desired targets.The solid curve represents population growth corresponding to the target set T1 and the dotted curve corresponds to the set T2 .It is clear from this result that population increases much more rapidly during the first 5-6 years.This is because the target is far from the initial population and the error is large which encourages maximum admissible immigration rate for reaching the lower boundary as quickly as possible.The speed of approach is dependent on the maximum admissible immigration rate uM .Once the total population exceeds the lower boundary, the growth slows down and the optimum immigration policy seems to maintain the population in the neighborhood of the target set.In contrast, if the population exceeds the upper boundary the optimum policy seems to pull it down by cutting down the immigration rate.
Intuitively these results suggest that if the current population is far below the desired target set, the optimum immigration rate should be set at its maximum permissible level so that the required manpower for economic activities is met as fast as possible.
Comparing the population growth (Fig. 12) corresponding to optimum policies with those of (Fig. 8 and Fig. 9) corresponding to the official immigration policy, it is clear that the total population can be well regulated and steered to any specified target.
(B) Variable Target  If it is required to reach the desired population level smoothly (over several years), the target set can be modified by use of a pair of smooth curves describing the upper and the lower boundaries starting around the initial population see Fig. 14.For illustration, we chose the variable target set given by The results are shown in Fig. 13-14.Fig. 13 shows the optimal control (immigration) policy.It is clear from the graph that it is smooth (not bang-bang as in the case of fixed target set) and the corresponding population grows smoothly and remains confined in the target set.
We conclude from these results that control theory provides a promising tool for determining the optimum immigration policy seeking specified targets.In fact one can add to the cost function as many factors as one desires, including humanitarian factors, to reflect the concerns of the society and use the optimization methodology proposed here to determine the optimum policy.This technique can be used by the Department of Immigration as an intelligent tool for determining the optimum immigration policy.
We must mention that the numerical results presented here are applicable only for the years 1976-2002.For current and future years one must use the corresponding data for numerical simulation and optimization and derive the optimum policies.

Immigration Policy Versus Job Creation Rate
Although many factors may affect the immigration policy, the unemployment rate is one of the most important one (Siklos and Marr, 1998:127-147;Denton, Feaver and Spencer,1997;Veugelers and Klassen, 1994: 351-370).The question is what should be the intake rate (immigration rate) so as to satisfy the manpower demand and at the same time keep the unemployment rate low.In addition to humanitarian factors, it is reasonable to tie the immigration rate with availability of jobs and job creation rate.Otherwise one can expect many social and political problems.16) of the previous section does not include the unemployment factor.To introduce this factor, we note that the growth or decline of unemployment rate is proportional to the growth or decline of population of age group G2 (the employable population) and the job growth rate.For this purpose we introduce a fourth state variable x4 , the unemployed population, and modify the model ( 16) to the following one, and the gradient vector is given by For numerical results, using the state equation ( 26) and the co-state equation ( 31) along with the gradient (34) as defined above, we use the same algorithm as indicated in section 4.1.

Unemployment Rate
According to Stat.Canada definition, the unemployment rate is given by the ratio of the number of unemployed people in labor force to the number of people in the labor force.Here, we only consider labor force in the age group G2.
For composition of labor force, a certain percentage of the population in age group G2 is eliminated.These are people in the age group G2 who are not actively seeking for jobs (such as students, handicaps, terminally sick, etc.).According to Stat.Canada definition, the labor force is given by some percentage of this population.That is, The factor p is a function of time as shown in Figure 15.The unemployment rate is then expressed by:

Numerical Results
Data Used For numerical results, the data used are as given below: Figure 16 shows the optimum immigration rate (or intake rate) of population of age group G2.Again this is approximately 500,000 per annum during the early period of the planning horizon.The curves 16A-B represent the numerical results corresponding to the Cases A-B respectively.The optimum immigration policies corresponding to the cases A and B are very close and hence the total population of the two cases, which is shown in Figure 18, coincides.But the optimum job creation rates are significantly different.Since job creation rate in case B is higher compared to that in case A, the unemployment rate in case B is lower.
Again it is clear that for fixed target set (time invariant), the optimal control policy is nearly bang-bang.As a result of this control policy, the population rapidly increases initially and then slows down preventing large deviation from the upper boundary of the target set.Again, we expect control policies to be smooth if variable target set T(t) is chosen as in subsection 5.2.

Conclusion
In this paper we have demonstrated that by use of modern Systems and Optimal Control theory, it is possible to formulate optimum immigration and job creation strategies while maintaining population level close to certain pre-specified target sets.
We have constructed in section 2, a simplified dynamic population model using the models proposed in (Ahmed and Rahim, 2001: 325-358).Using the basic data (birth, death and transition rates etc.) from Stat.Canada, we found the numerical results (population) based on our model in close agreement with the actual population.This was presented in section 3.
In case of non availability of the basic parameters, one must use available population statistics to determine the unknown parameters.This is known as

283
identification, and we have demonstrated in section 4 that optimal control theory can be used to determine these parameters.
Based on the model constructed above, in section 5, we have formulated a control problem with the objective of reaching a specified population target set (fixed as well as variable) by use of immigration rate as the control variable.Optimal control theory is used to determine the optimum immigration policy as illustrated by numerical results.
In section 6, the population model is augmented by including a fourth equation describing the dynamics of unemployment rate, including job creation rate as another control variable.Following our methodology, optimum immigration and job creation policies were determined.Results are illustrated by numerical simulation and they are found to be very encouraging.
Using actual field data along with the desired objective functional, reflecting important social concerns, and following the methodology presented here one can develop optimum immigration and job creation policies.
rates from G1 to G2 and from G2 to G3, respectively.The parameter 3 denotes the birth rate in group G1, and the three parameters 4 , 5 and 6 denote the death rates in groups G1,G2, and G3 respectively.The initial population (condition) of the three age groups in the year 1972, is given by: cost function as being the identification error, denotes the population vector obtained by solving the model equation (7) corresponding to any arbitrary choice of the parameter .population obtained from census data for the period [0, T].
are obtained by using the mean of actual values reported by Statistics Canada.This means that 31percent of the number of adults immigrated are children.The range of the control variable is defined by the limit means that immigration rate must not exceed 0.05 percent of the adult population G2.This limit can be fixed by the immigration department on the basis of other national concerns.The initial population of the three age groups in the year 1976, is given by: xm , xM denote the lower and upper boundaries of the target population [ xm , xM ], 1 and 2 are the weights assigned to deviations from the boundaries of the target set.They are chosen (A) Initial Condition for the State EquationThe initial condition for the state equation is given by: For better illustration and comparison of results we use two different sets of control constraints (A) and (B).The range of control variables u u 2 1 , are described in Table1.The upper limit of the job creation rate in case (B) is larger than that of case (A).(C) Percentage of Labor force (Fig.15) Percentage of labor force, which is obtained from Stat.Canada, is shown in Figure 15.Figures 16-19. 262) For notational convenience, we use A to denote the sum of terms as shown