# A Mathematical Overview of the Monogamous Marriage in a Multiregion Framework: Modelling and Control

Authors: Mustapha Lhous; Omar Zakary; Mostafa Rachik
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

The main objective of this paper is to develop a new mathematical model to study, analyze, and control the family status in several regions and to discuss the impact of the connectivity of regions and the mobility of residents on the marital status of the family, by adopting a multiregion discrete-time model. The modelling and the control process of the system that describes the case of monogamous marriages in a multiregion framework are considered. Two combined control strategies are proposed, which allow reducing the virgin and divorced individuals and increasing the number of married individuals in a specific region. The first control is considered as the impact of public awareness campaigns to educate virgin men and women about the benefits of marriage for the individual and the society; the second control characterizes the legal procedures, administrative complications, and the heavy financial and social consequences of divorces. The optimal control theory is applied to characterize such optimal strategies and determined numerically using a progressive-regressive discrete scheme to discuss the obtained results.


---

## Body

## 1. Introduction

Mathematical models are useful tools for understanding the functioning of natural systems and for predicting their evolution. Among these models are those that study the dynamics of populations and ecosystems. Many researchers have studied models of population dynamics: prey and predator dynamics [1–4], epidemic dynamics in a population [5–8], molecular systems [9–11], and so forth.Civil status is the situation of the person in the family and/or society. A person’s marital status is positioned in one of four categories: virgin, married, divorced, or widowed (VMDW). Several authors have discussed the social functioning, marital status, family stability, and social control of health behavior [12–18]. In [19], the authors have given a discrete mathematical model that describes the marital status of the monogamous marriage case and they also studied the optimal control that reduces the number of unmarried and divorced people and increases the number of married people. In this paper, we generalize the results established in [19] to the case of a multiregion system. We studied the effect of the population travel to different regions for the marriage, divorce, and widow.Many studies on civil status have given statistics on the family situation of different regions (cities, rural regions, industrial regions, tourist regions, etc.), according to the age group of the population. For example in Morocco, the High Commissioner for Planning (HCP, public institution) [14] gives different statistics according to the age groups of the population; there are 76.2% of single people for the age group of 20 to 24 and 54% for the group of 25 to 29 years old. The percentage of divorce is 30% for couples with less than 5 years of marriage and 3% for married couples who have exceeded 20 years of marriages. 41.4% of married men and 41.8% of married women are identified by the HCP in the urban areas. The percentage of divorce is 0.8% for men and 3.1% for women. The widowed men are 0.6% of the population and women are 7.3% in urban areas. For the rural area, there are 40% of married men and 42.5% of married women; the percentage of divorce is 0.4% for married men and 1.4% for married women. Widowed men are represented by 0.6% of the population and women by 6.8% in rural areas.In recent years, many attempts have been made to develop control strategies for different systems [5, 20–27]. There are a number of different methods for calculating the optimal control for a specific mathematical model. For example, Pontryagin et al.’s maximum principle [28] allows the characterization of the optimal control for an ordinary equation model system with a given constraint.In [29], the authors have described a new modelling approach based on multiregion discrete-time SIR models aiming to describe the spatial-temporal evolution of epidemics that emerge in different geographical regions and to show the influence of one region on another region via infection travel.In this paper and inspired by [29], we investigate an approach that determines an optimal control relative to a discrete VMDW model in a multiregion framework, which defines the evolution of the marital status of the marriage in a population, enabling decision-makers to develop very useful control strategies to reduce the virgin individuals and to increase the number of married individuals in a specific region.The first control can be considered as public awareness campaigns showing to individuals the benefits of marriage on the psychic and social stability on persons and the society, or cultural entertained events to give people the chance to meet and to know and to allow themselves to get married. The second control is determined for the persons who have initiated the divorce proceedings; this control is considered as a long and costly legal procedure.The optimal control problem was the subject of an optimization criterion represented by the minimization of an objective function. The optimality system is solved based on an iterative discrete schema that converges following an appropriate test similar to the one related to the Forward-Backward Sweep Method (FBSM).The paper is organized as follows: in Section2, the model VMDW is described for a multiregion discrete model. In Section 3, we give some results concerning the existence of optimal control and we use Pontryagin’s maximum principle to investigate the analysis of control strategies and to determine the necessary condition for optimal control. Numerical simulations are given in Section 4. Finally, we conclude the paper in Section 5.


## 2. Mathematical Modelling

We consider a discrete-time model VMDW of the marital status of the family dynamic within a domain of interestΩ which represents a country, a city, a town, or a small domain. We assume that there are p geographical regions (domains) of the domain studied Ω. Let Ω=∪j=1pΩj , and let NiΩk be the population of domain Ωk at time i ; that is, the number of individuals who are residents in Ωk.This model classifies the marital status of the family dynamics of a population into eight compartments in each regionΩk: virgin men VMΩk, virgin women VWΩk, married men MMΩk, married women MWΩk, divorced men DMΩk, divorced women DWΩk, widowed men WMΩk, and widowed women WWΩk.The unit of timei can correspond to days, months, or years; it depends on the frequency of the survey and demographic studies as needed. However demographic statistics are generally done annually so the units i, i+1,… can be considered as years. The following system describes a discrete model of the marital status of the monogamous marriage case of a region Ωk:(1)Vi+1MΩk=Λk1+ViMΩk−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩjViMΩk−dkViMΩk,(2)Vi+1WΩk=Λk2+ViWΩk−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩkViWΩk−dkViWΩk,(3)Di+1MΩk=DiMΩk−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩjDiMΩk+λkMiMΩk−dkDiMΩk,(4)Di+1WΩk=DiWΩk−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩkDiWΩk+∑j=1pλjkMiWΩj−dkDiWΩk,(5)Wi+1MΩk=WiMΩk−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩjWiMΩk+ωkMiWΩk−dkWiMΩk,(6)Wi+1WΩk=WiWΩk−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩkWiWΩk+∑j=1pρjkMiMΩj−dkWiWΩk,(7)Mi+1MΩk=MiMΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩjViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ρkMiMΩk,(8)Mi+1WΩk=MiWΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ωkMiWΩk,where V0MΩk, V0WΩk, M0MΩk, M0WΩk, D0MΩk, D0WΩk, W0MΩk, and W0WΩk are the given initial state in the region Ωk. In equations (1)–(8), all parameters are nonnegative and defined in Table 1.Table 1
The description of parameters used for the definition of discrete-time systems (1)–(8).
ParameterDescriptionΛk1Natality of virgin men in regionΩkΛk2Natality of virgin women in regionΩkdkNatural mortality rate in regionΩkαkjMarriage rate of virgin men of regionΩk to virgin women of region ΩjβkjMarriage rate of divorced men of regionΩk to virgin women of region ΩjγkjMarriage rate of virgin men of regionΩk to divorced women of region ΩjδkjMarriage rate of virgin men of regionΩk to widowed women of region ΩjηkjMarriage rate of widowed men of regionΩk to virgin women of region ΩjμkjMarriage rate of divorced men of regionΩk to divorced women of region ΩjθkjMarriage rate of divorced men of regionΩk to widowed women of region ΩjνkjMarriage rate of widowed men of regionΩk to divorced women of region ΩjσkjMarriage rate of widowed men of regionΩk to widowed women of region ΩjλkDivorce rate of married men of regionΩkλjkDivorce rate of women of regionΩk which return to the region ΩjρkWidow rate of married women of regionΩkωkWidow rate of married men of regionΩkFor equation (1) of the model, virgin men of a region Ωk can contact a virgin, divorced, or widowed woman of another region Ωj with αkj, γkj, and δkj rates, respectively, and this contact can result in a marriage. The married couple can stay in the region Ωk or settle in the conjoin region Ωj. Thus, the number of virgin men decreases and the number of virgins at the instant i is substituted for the number αkjViMΩkViWΩj/NiΩk+NiΩj+γkjViMΩkDiWΩj/NiΩk+NiΩj+δkjViMΩkWiWΩj/NiΩk+NiΩj and a portion of this number αkjkViMΩkViWΩj/NiΩk+NiΩj+γkjkViMΩkDiWΩj/NiΩk+NiΩj+δkjkViMΩkWiWΩj/NiΩk+NiΩj that represents the part of married couples who decides to settle in Ωk is added at the time i+1 to the number of men marrying in region Ωk and the other portion αkjjViMΩkViWΩj/NiΩk+NiΩj+γkjjViMΩkDiWΩj/NiΩk+NiΩj+δkjjViMΩkWiWΩj/NiΩk+NiΩj is added to the number of men marrying in region Ωj. Then, we have αkj=αkjk+αkjj, γkj=γkjk+γkjj and δkj=δkjk+δkjj.Similarly in equation (2), the number of virgin women in region Ωk decreases at the instant i+1 by substituting the number of virgin women at the instant i, the number αjkViMΩjViWΩk/NiΩj+NiΩk+βjkDiMΩjViWΩk/NiΩj+NiΩk+ηjkWiMΩjViWΩk/NiΩj+NiΩk which represents the number of married women after contact with a virgin, divorced, or widowed men in region Ωj with αjk, βjk, and ηjk levels, respectively. Also we have αjk=αjkj+αjkk, βjk=βjkj+βjkk and ηjk=ηjkj+ηjkk.In equation (3), a divorced man in region Ωk can contact a virgin, divorced, or widowed woman of region Ωj with βkj, μkj, and θkj rates, respectively, so this contact can result in marriage. Then the number βkjDiMΩkViWΩj/NiΩk+NiΩj+μkjDiMΩkDiWΩj/NiΩk+NiΩj+θkjDiMΩkWiWΩj/NiΩk+NiΩj is decreased by the number of divorced men at the time i+1 and added to the number of married men and we have μkj=μkjk+μkjj and θkj=θkjk+θkjj. In the model we propose here, it was considered that a divorced man remains in his region and the divorced woman returns to the region of her parents; this is the case for the majority of regions whether they are conservative or not. And so the number λkMiMΩk is added to the number of divorced men and the number λjkMiWΩj of divorced women is added to the number of divorced women of the region Ωk with λjk being the rate that a woman divorces a man from the region Ωk and returns to the region Ωj. Consequently, the total number of ∑j=1pλjkMiWΩj is added to the number of divorced women in the region Ωk with the relation λk=∑j=1pλkj. The same principle applies to equations (5) and (6) with ρk=∑j=1pρkj.For equation (7), the number of married men in the region Ωk increases at the instant i+1 by the number of virgin, divorced, and widowed men who are married by contacting virgin, divorced, or widowed women of the region Ωj and decreases the natural mortality with a ρk and divorce rate with a λk rate. The same principle can be applied to equation (8).In equations (7) and (8) from the model, the natural mortality of married men and women was considered. The mortality rate for married men is ρk and the mortality rate for married women is ωk. These rates appear in equations (5) and (6) which correspond to widowed men and women, respectively, so the number of mortalities which is ρkMiMΩk married men and the number ∑j=1pρjkMiMΩj are added to the number of widowed women and similarly the number ωkMiWΩk is added to the number of widowed men.


## 3. Methods and Results

### 3.1. An Optimal Control Approach

An optimal control approach has been applied to models (1)–(8) to reduce the virgin and divorced individuals and increase the number of married individuals along the control strategy period. For this, we introduce a control variable ui,vi that characterizes the benefits of an awareness campaign to educate virgin men and women about the benefits of marriage for the individual and the society, especially the legal procedures, administrative complications, and the heavy financial and social consequences of divorces, respectively, in the abovementioned models (1)–(8). Then, the model is given by the following equations:(9)Vi+1MΩk=Λk1+ViMΩk−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩjViMΩk−dkViMΩk−uiViMΩk,(10)Vi+1WΩk=Λk2+ViWΩk−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩkViWΩk−dkViWΩk−uiViWΩk,(11)Di+1MΩk=DiMΩk−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩjDiMΩk+λkMiMΩk−dkDiMΩk−viDiMΩk,(12)Di+1WΩk=DiWΩk−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩkDiWΩk+∑j=1pλjkMiWΩj−dkDiWΩk−viDiWΩk,(13)Wi+1MΩk=WiMΩk−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩjWiMΩk+ωkMiWΩk−dkWiMΩk,(14)Wi+1WΩk=WiWΩk−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩkWiWΩk+∑j=1pρjkMiMΩj−dkWiWΩk,(15)Mi+1MΩk=MiMΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ρkMiMΩk+uiViMΩk+viDiMΩk,(16)Mi+1WΩk=MiWΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ωkMiWΩk+uiViWΩk+viDiWΩk,


### 3.2. Characterization of the Optimal Control

For an initial stateV0MΩk,V0WΩk,M0MΩk,M0WΩk,D0MΩk,D0WΩk, W0MΩk,W0WΩk, we consider an optimization criterion defined by the following objective function:(17)Ju,v=∑i=0NA1ViWΩk+A2DiWΩk−A3MiWΩk+∑i=0N−1τ12ui2+τ22vi2,subject to systems (9)–(16). Here, A1, A2, and A3 are positive constants to keep a balance in the size of ViWΩk, DiWΩk, and MiWΩk, respectively. In the objective functional, τ1 and τ2 are the positive weight parameters which are associated with the controls ui and vi.In other words, we seek the optimal controlsu∗,v∗ such that(18)Ju∗,v∗=minJu,vu,v∈Uad,where Uad is the set of admissible controls defined by(19)Uad=u,vumin≤ui≤umax,vmin≤vi≤vmax,i∈0,…,N−1,where umin,umax,vmin,vmax∈0,14.The sufficient condition for existence of an optimal controlu∗,v∗ for problem (18) follows from standard results of [21]. In order to find an optimal solution, first we find the Hamiltonian for the optimal control problem (18). In fact, the Hamiltonian H of the optimal problem is given by(20)H=A1ViWΩk+A2DiWΩk−A3MiWΩk+τ12ui2+τ22vi2+ζ1,i+1kΛk1+ViMΩk−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩjViMΩk−dkViMΩk−uiViMΩk+ζ2,i+1kΛk2+ViWΩk−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩkViWΩk−dkViWΩk−uiViWΩk+ζ3,i+1kDiMΩk−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩjDiMΩk+λkMiMΩk−dkDiMΩk−viDiMΩk+ζ4,i+1kDiWΩk−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩkDiWΩk+∑j=1pλjkMiWΩj−dkDiWΩk−viDiWΩk+ζ5,i+1kWiMΩk−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩjWiMΩk+ωkMiWΩk−dkWiMΩk+ζ6,i+1kWiWΩk−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩkWiWΩk+∑j=1pρjkMiMΩj−dkWiWΩk+ζ7,i+1kMiMΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ρkMiMΩk+uiViMΩk+viDiMΩk+ζ8,i+1kMiWΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ωkMiWΩk+uiViWΩk+viDiWΩk,where ζ1k,ζ2k,…,ζ8k are the adjoint functions to be determined suitably.At the same time by using Pontryagin et al.’s maximum principle [28], we derive necessary conditions for our optimal control. We obtain the following theorem.Theorem 1.
LetVM∗Ωk, VW∗Ωk, MM∗Ωk, MW∗Ωk, DM∗Ωk, DW∗Ωk, WM∗Ωk, and WW∗Ωk be optimal state solutions with associated optimal control u∗,v∗ for the optimal control problem (18). Then, there exist adjoint variables ζ1k,ζ2k,…,ζ8k that satisfy(21)Δζ1,ik=−ζ1,i+1k1−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩj−dk−ui+ζ2,i+1kαkkViWΩk2NiΩk+ζ4,i+1kγkkDiWΩk2NiΩk+ζ6,i+1kδkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkkViWΩk+γkkkDiWΩk+δkkkWiWΩk2NiΩk−ζ7,i+1kui,(22)Δζ2,ik=−A1−ζ2,i+1k1−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩk−dk−ui+ζ1,i+1kαkkViMΩk2NiΩk+ζ3,i+1kβkkDiMΩk2NiΩk+ζ5,i+1kηkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩk+αkkkViMΩk+βkkkDiMΩk+ηkkkWiMΩk2NiΩk−ζ8,i+1kui,(23)Δζ3,ik=−ζ3,i+1k1−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩj−dk−vi+ζ2,i+1kβkkViWΩk2NiΩk+ζ4,i+1kμkkDiWΩk2NiΩk+ζ6,i+1kθkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+βkkkViWΩk+μkkkDiWΩk+θkkkWiWΩk2NiΩk−ζ7,i+1kvi,(24)Δζ4,ik=−A2−ζ4,i+1k1−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩk−dk−vi+ζ1,i+1kγkkViMΩk2NiΩk+ζ3,i+1kμkkDiMΩk2NiΩk+ζ5,i+1kνkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩk+γkkkViMΩk+μkkkDiMΩk+νkkkWiMΩk2NiΩk−ζ8,i+1kvi,(25)Δζ5,i=−ζ5,i+1k1−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩj−dk+ζ2,i+1kηkkViWΩk2NiΩk+ζ4,i+1kνkkDiWΩk2NiΩk+ζ6,i+1kσkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+ηkkkViWΩk+νkkkDiWΩk+σkkkWiWΩk2NiΩk,(26)Δζ6,ik=−ζ6,i+1k1−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩk−dk+ζ1,i+1kδkkViMΩk2NiΩk+ζ3,i+1kθkkDiMΩk2NiΩk+ζ5,i+1kσkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩk+δkkkViMΩk+θkkkDiMΩk+σkkkWiMΩk2NiΩk,(27)Δζ7,ik=−ζ3,i+1kλk−ζ6,i+1kρkk−ζ7,i+1k1−λk−ρk,(28)Δζ8,ik=A3−ζ4,i+1kλkk−ζ5,i+1kωk−ζ8,i+1k1−λk−ωk,with transversality conditions(29)ζ1,N=0,ζ2,N=A1,ζ3,N=0,ζ4,N=A2,ζ5,N=0,ζ6,N=0,ζ7,N=0,ζ8,N=−A3.
Furthermore, the optimal controlui∗,vi∗ is given by(30)ui∗=minmaxViMΩkζ1,i+1−ζ7,i+1+ViWΩkζ2,i+1−ζ8,i+1τ1,umin,umax.(31)vi∗=minmaxDiMΩkζ3,i+1−ζ7,i+1+DiWΩkζ4,i+1−ζ8,i+1τ2,vmin,vmax,for i=0,…,N−1.Proof.
Using Pontryagin et al.’s maximum principle [28] and setting VM∗Ωk, VW∗Ωk, MM∗Ωk, MW∗Ωk, DM∗Ωk, DW∗Ωk, WM∗Ωk, WW∗Ωk, and u∗,v∗, we obtain the following adjoint equations:(32)Δζ1,ik=−∂H∂ViMΩk=−ζ1,i+11−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩj−dk−ui−ζ2,i+1kαkk2NiΩkViWΩk−ζ4,i+1kγkk2NiΩkDiWΩk−ζ6,i+1kδkk2NiΩkWiWΩk+ζ7,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkk2NiΩkViWΩk+γkkk2NiΩkDiWΩk+δkkk2NiΩkWiWΩk+ui+ζ8,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkk2NiΩkViWΩk+γkkk2NiΩkDiWΩk+δkkk2NiΩkWiWΩk=−ζ1,i+1k1−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩj−dk−ui+ζ2,i+1kαkkViWΩk2NiΩk+ζ4,i+1kγkkDiWΩk2NiΩk+ζ6,i+1kδkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkkViWΩk+γkkkDiWΩk+δkkkWiWΩk2NiΩk−ζ7,i+1kui,(33)Δζ2,ik=−∂H∂ViWΩk=−A1−ζ1,i+1kαkk2NiΩkViMΩkζ2,i+11−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩk+NiΩj−dk−ui−ζ3,i+1kβkk2NiΩkDiMΩk−ζ5,i+1kηkk2NiΩkWiMΩk+ζ7,i+1kαkkk2NiΩkViMΩk+∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩk+NiΩj+βkkk2NiΩkDiMΩk+ηkkk2NiΩkWiMΩk+ζ8,i+1kαkkk2NiΩkViMΩk+∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩk+NiΩj+βkkk2NiΩkDiMΩk+ηkkk2NiΩkWiMΩk+ui=−A1−ζ2,i+1k1−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩk−dk−ui+ζ1,i+1kαkkViMΩk2NiΩk+ζ3,i+1kβkkDiMΩk2NiΩk+ζ5,i+1kηkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩk+αkkkViMΩk+βkkkDiMΩk+ηkkkWiMΩk2NiΩk−ζ8,i+1kui,(34)Δζ3,ik=−∂H∂DiMΩk=−−ζ2,i+1kβkk2NiΩkViWΩk+ζ3,i+11−∑j=1pβkjViWΩj+μjkDiWΩj+θkjWiWΩjNiΩk+NiΩj−dk−vi−ζ4,i+1kμkk2NiΩkDiWΩk−ζ6,i+1kθkk2NiΩkWiWΩk+ζ7,i+1kβkkk2NiΩkViWΩk+∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+μkkk2NiΩkDiWΩk+θkkk2NiΩkWiWΩk+vi+ζ8,i+1kβkkk2NiΩkViWΩk+∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+μkkk2NiΩkDiWΩk+θkkk2NiΩkWiWΩk=−ζ3,i+1k1−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩj−dk−vi+ζ2,i+1kβkkViWΩk2NiΩk+ζ4,i+1kμkkDiWΩk2NiΩk+ζ6,i+1kθkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+βkkkViWΩk+μkkkDiWΩk+θkkkWiWΩk2NiΩk−ζ7,i+1kvi,(35)Δζ4,ik=−∂H∂DiWΩk=−A2−ζ1,i+1kγkk2NiΩkViMΩk−ζ3,i+1kμkk2NiΩkDiMΩk+ζ4,i+11−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩk+NiΩj−dk−vi−ζ5,i+1kνkk2NiΩkWiMΩk+ζ7,i+1kγkk2NiΩkViMΩk+μkkk2NiΩkDiMΩk+∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩk+NiΩj+νkkk2NiΩkWiMΩk+ζ8,i+1kγkk2NiΩkViMΩk+μkkk2NiΩkDiMΩk+∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩk+NiΩj+νkkk2NiΩkWiMΩk+vi=−A2−ζ4,i+1k1−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩk−dk−vi+ζ1,i+1kγkkViMΩk2NiΩk+ζ3,i+1kμkkDiMΩk2NiΩk+ζ5,i+1kνkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩk+γkkkViMΩk+μkkkDiMΩk+νkkkWiMΩk2NiΩk−ζ8,i+1kvi,(36)Δζ5,ik=−∂H∂WiMΩk=−−ζ2,i+1kηkk2NiΩkViWΩk−ζ4,i+1kνkk2NiΩkDiWΩk+ζ5,i+11−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩj−dk−ζ6,i+1kσkk2NiΩkWiWΩk+ζ7,i+1kηkk2NiΩkViWΩk+νkkk2NiΩkDiWΩk+∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+σkkk2NiΩkWiWΩk+ζ8,i+1kηkk2NiΩkViWΩk+νkkk2NiΩkDiWΩk+∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+σkkk2NiΩkWiWΩk=−ζ5,i+1k1−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩj−dk+ζ2,i+1kηkkViWΩk2NiΩk+ζ4,i+1kνkkDiWΩk2NiΩk+ζ6,i+1kσkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+ηkkkViWΩk+νkkkDiWΩk+σkkkWiWΩk2NiΩk,(37)Δζ6,ik=−∂H∂WiWΩk=−−ζ1,i+1kδkk2NiΩkViMΩk−ζ3,i+1kθkk2NiΩkDiMΩk−ζ5,i+1kσkk2NiΩkWiMΩk+ζ6,i+11−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩk+NiΩj−dk+ζ7,i+1kδkk2NiΩkViMΩk+θkkk2NiΩkDiMΩk+σkkk2NiΩkWiMΩk+∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩk+NiΩj+ζ8,i+1kδkk2NiΩkViMΩk+θkkk2NiΩkDiMΩk+σkkk2NiΩkWiMΩk+∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩk+NiΩj=−ζ6,i+1k1−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩk−dk+ζ1,i+1kδkkViMΩk2NiΩk+ζ3,i+1kθkkDiMΩk2NiΩk+ζ5,i+1kσkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩk+δkkkViMΩk+θkkkDiMΩk+σkkkWiMΩk2NiΩk,(38)Δζ7,ik=−∂H∂MiMΩk=−ζ3,i+1kλk+ζ6,i+1kρkk+ζ7,i+1k1−λk+ρk,(39)Δζ8,ik=−∂H∂MiWΩk=−−A3+ζ4,i+1kλk+ζ5,i+1kωk+ζ8,i+1k1−λk+ωk,with transversality conditions:(40)ζ1,N=0,ζ2,N=A1,ζ3,N=0,ζ4,N=A2,ζ5,N=0,ζ6,N=0,ζ7,N=0ζ8,N=−A3.
To obtain the optimality conditions, we take the variation with respect to controlui∗,vi∗ and set it equal to zero:(41)∂H∂ui=τ1ui−ζ1,i+1kViMΩk−ζ2,i+1kViWΩk+ζ7,i+1kViMΩk+ζ8,i+1kViWΩk=0,∂H∂vi=τ2ui−ζ3,i+1kDiMΩk−ζ4,i+1kDiWΩk+ζ7,i+1kDiMΩk+ζ8,i+1kDiWΩk=0.
Then, we obtain the optimal control:(42)ui∗=ViMΩkζ1,i+1k−ζ7,i+1k+ViWΩkζ2,i+1k−ζ8,i+1kτ1;vi∗=DiMΩkζ3,i+1k−ζ7,i+1k+DiWΩkζ4,i+1k−ζ8,i+1kτ2.
By the bounds inUad, it is easy to obtain ui∗,vi∗ in the following form:(43)ui∗=minmaxViMΩkζ1,i+1k−ζ7,i+1k+ViWΩkζ2,i+1k−ζ8,i+1kτ1,umin,umax,vi∗=minmaxDiMΩkζ3,i+1k−ζ7,i+1k+DiWΩkζ4,i+1k−ζ8,i+1kτ2,umin,umax,for i=0,…,N−1.


## 3.1. An Optimal Control Approach

An optimal control approach has been applied to models (1)–(8) to reduce the virgin and divorced individuals and increase the number of married individuals along the control strategy period. For this, we introduce a control variable ui,vi that characterizes the benefits of an awareness campaign to educate virgin men and women about the benefits of marriage for the individual and the society, especially the legal procedures, administrative complications, and the heavy financial and social consequences of divorces, respectively, in the abovementioned models (1)–(8). Then, the model is given by the following equations:(9)Vi+1MΩk=Λk1+ViMΩk−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩjViMΩk−dkViMΩk−uiViMΩk,(10)Vi+1WΩk=Λk2+ViWΩk−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩkViWΩk−dkViWΩk−uiViWΩk,(11)Di+1MΩk=DiMΩk−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩjDiMΩk+λkMiMΩk−dkDiMΩk−viDiMΩk,(12)Di+1WΩk=DiWΩk−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩkDiWΩk+∑j=1pλjkMiWΩj−dkDiWΩk−viDiWΩk,(13)Wi+1MΩk=WiMΩk−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩjWiMΩk+ωkMiWΩk−dkWiMΩk,(14)Wi+1WΩk=WiWΩk−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩkWiWΩk+∑j=1pρjkMiMΩj−dkWiWΩk,(15)Mi+1MΩk=MiMΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ρkMiMΩk+uiViMΩk+viDiMΩk,(16)Mi+1WΩk=MiWΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ωkMiWΩk+uiViWΩk+viDiWΩk,


## 3.2. Characterization of the Optimal Control

For an initial stateV0MΩk,V0WΩk,M0MΩk,M0WΩk,D0MΩk,D0WΩk, W0MΩk,W0WΩk, we consider an optimization criterion defined by the following objective function:(17)Ju,v=∑i=0NA1ViWΩk+A2DiWΩk−A3MiWΩk+∑i=0N−1τ12ui2+τ22vi2,subject to systems (9)–(16). Here, A1, A2, and A3 are positive constants to keep a balance in the size of ViWΩk, DiWΩk, and MiWΩk, respectively. In the objective functional, τ1 and τ2 are the positive weight parameters which are associated with the controls ui and vi.In other words, we seek the optimal controlsu∗,v∗ such that(18)Ju∗,v∗=minJu,vu,v∈Uad,where Uad is the set of admissible controls defined by(19)Uad=u,vumin≤ui≤umax,vmin≤vi≤vmax,i∈0,…,N−1,where umin,umax,vmin,vmax∈0,14.The sufficient condition for existence of an optimal controlu∗,v∗ for problem (18) follows from standard results of [21]. In order to find an optimal solution, first we find the Hamiltonian for the optimal control problem (18). In fact, the Hamiltonian H of the optimal problem is given by(20)H=A1ViWΩk+A2DiWΩk−A3MiWΩk+τ12ui2+τ22vi2+ζ1,i+1kΛk1+ViMΩk−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩjViMΩk−dkViMΩk−uiViMΩk+ζ2,i+1kΛk2+ViWΩk−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩkViWΩk−dkViWΩk−uiViWΩk+ζ3,i+1kDiMΩk−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩjDiMΩk+λkMiMΩk−dkDiMΩk−viDiMΩk+ζ4,i+1kDiWΩk−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩkDiWΩk+∑j=1pλjkMiWΩj−dkDiWΩk−viDiWΩk+ζ5,i+1kWiMΩk−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩjWiMΩk+ωkMiWΩk−dkWiMΩk+ζ6,i+1kWiWΩk−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩkWiWΩk+∑j=1pρjkMiMΩj−dkWiWΩk+ζ7,i+1kMiMΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ρkMiMΩk+uiViMΩk+viDiMΩk+ζ8,i+1kMiWΩk+∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩjViMΩk+αjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩkViWΩk+βkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩjDiMΩk+γjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩkDiWΩk+ηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩjWiMΩk+δjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩkWiWΩk−λk+ωkMiWΩk+uiViWΩk+viDiWΩk,where ζ1k,ζ2k,…,ζ8k are the adjoint functions to be determined suitably.At the same time by using Pontryagin et al.’s maximum principle [28], we derive necessary conditions for our optimal control. We obtain the following theorem.Theorem 1.
LetVM∗Ωk, VW∗Ωk, MM∗Ωk, MW∗Ωk, DM∗Ωk, DW∗Ωk, WM∗Ωk, and WW∗Ωk be optimal state solutions with associated optimal control u∗,v∗ for the optimal control problem (18). Then, there exist adjoint variables ζ1k,ζ2k,…,ζ8k that satisfy(21)Δζ1,ik=−ζ1,i+1k1−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩj−dk−ui+ζ2,i+1kαkkViWΩk2NiΩk+ζ4,i+1kγkkDiWΩk2NiΩk+ζ6,i+1kδkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkkViWΩk+γkkkDiWΩk+δkkkWiWΩk2NiΩk−ζ7,i+1kui,(22)Δζ2,ik=−A1−ζ2,i+1k1−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩk−dk−ui+ζ1,i+1kαkkViMΩk2NiΩk+ζ3,i+1kβkkDiMΩk2NiΩk+ζ5,i+1kηkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩk+αkkkViMΩk+βkkkDiMΩk+ηkkkWiMΩk2NiΩk−ζ8,i+1kui,(23)Δζ3,ik=−ζ3,i+1k1−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩj−dk−vi+ζ2,i+1kβkkViWΩk2NiΩk+ζ4,i+1kμkkDiWΩk2NiΩk+ζ6,i+1kθkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+βkkkViWΩk+μkkkDiWΩk+θkkkWiWΩk2NiΩk−ζ7,i+1kvi,(24)Δζ4,ik=−A2−ζ4,i+1k1−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩk−dk−vi+ζ1,i+1kγkkViMΩk2NiΩk+ζ3,i+1kμkkDiMΩk2NiΩk+ζ5,i+1kνkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩk+γkkkViMΩk+μkkkDiMΩk+νkkkWiMΩk2NiΩk−ζ8,i+1kvi,(25)Δζ5,i=−ζ5,i+1k1−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩj−dk+ζ2,i+1kηkkViWΩk2NiΩk+ζ4,i+1kνkkDiWΩk2NiΩk+ζ6,i+1kσkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+ηkkkViWΩk+νkkkDiWΩk+σkkkWiWΩk2NiΩk,(26)Δζ6,ik=−ζ6,i+1k1−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩk−dk+ζ1,i+1kδkkViMΩk2NiΩk+ζ3,i+1kθkkDiMΩk2NiΩk+ζ5,i+1kσkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩk+δkkkViMΩk+θkkkDiMΩk+σkkkWiMΩk2NiΩk,(27)Δζ7,ik=−ζ3,i+1kλk−ζ6,i+1kρkk−ζ7,i+1k1−λk−ρk,(28)Δζ8,ik=A3−ζ4,i+1kλkk−ζ5,i+1kωk−ζ8,i+1k1−λk−ωk,with transversality conditions(29)ζ1,N=0,ζ2,N=A1,ζ3,N=0,ζ4,N=A2,ζ5,N=0,ζ6,N=0,ζ7,N=0,ζ8,N=−A3.
Furthermore, the optimal controlui∗,vi∗ is given by(30)ui∗=minmaxViMΩkζ1,i+1−ζ7,i+1+ViWΩkζ2,i+1−ζ8,i+1τ1,umin,umax.(31)vi∗=minmaxDiMΩkζ3,i+1−ζ7,i+1+DiWΩkζ4,i+1−ζ8,i+1τ2,vmin,vmax,for i=0,…,N−1.Proof.
Using Pontryagin et al.’s maximum principle [28] and setting VM∗Ωk, VW∗Ωk, MM∗Ωk, MW∗Ωk, DM∗Ωk, DW∗Ωk, WM∗Ωk, WW∗Ωk, and u∗,v∗, we obtain the following adjoint equations:(32)Δζ1,ik=−∂H∂ViMΩk=−ζ1,i+11−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩj−dk−ui−ζ2,i+1kαkk2NiΩkViWΩk−ζ4,i+1kγkk2NiΩkDiWΩk−ζ6,i+1kδkk2NiΩkWiWΩk+ζ7,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkk2NiΩkViWΩk+γkkk2NiΩkDiWΩk+δkkk2NiΩkWiWΩk+ui+ζ8,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkk2NiΩkViWΩk+γkkk2NiΩkDiWΩk+δkkk2NiΩkWiWΩk=−ζ1,i+1k1−∑j=1pαkjViWΩj+γkjDiWΩj+δkjWiWΩjNiΩk+NiΩj−dk−ui+ζ2,i+1kαkkViWΩk2NiΩk+ζ4,i+1kγkkDiWΩk2NiΩk+ζ6,i+1kδkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαkjkViWΩj+γkjkDiWΩj+δkjkWiWΩjNiΩk+NiΩj+αkkkViWΩk+γkkkDiWΩk+δkkkWiWΩk2NiΩk−ζ7,i+1kui,(33)Δζ2,ik=−∂H∂ViWΩk=−A1−ζ1,i+1kαkk2NiΩkViMΩkζ2,i+11−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩk+NiΩj−dk−ui−ζ3,i+1kβkk2NiΩkDiMΩk−ζ5,i+1kηkk2NiΩkWiMΩk+ζ7,i+1kαkkk2NiΩkViMΩk+∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩk+NiΩj+βkkk2NiΩkDiMΩk+ηkkk2NiΩkWiMΩk+ζ8,i+1kαkkk2NiΩkViMΩk+∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩk+NiΩj+βkkk2NiΩkDiMΩk+ηkkk2NiΩkWiMΩk+ui=−A1−ζ2,i+1k1−∑j=1pαjkViMΩj+βjkDiMΩj+ηjkWiMΩjNiΩj+NiΩk−dk−ui+ζ1,i+1kαkkViMΩk2NiΩk+ζ3,i+1kβkkDiMΩk2NiΩk+ζ5,i+1kηkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pαjkkViMΩj+βjkkDiMΩj+ηjkkWiMΩjNiΩj+NiΩk+αkkkViMΩk+βkkkDiMΩk+ηkkkWiMΩk2NiΩk−ζ8,i+1kui,(34)Δζ3,ik=−∂H∂DiMΩk=−−ζ2,i+1kβkk2NiΩkViWΩk+ζ3,i+11−∑j=1pβkjViWΩj+μjkDiWΩj+θkjWiWΩjNiΩk+NiΩj−dk−vi−ζ4,i+1kμkk2NiΩkDiWΩk−ζ6,i+1kθkk2NiΩkWiWΩk+ζ7,i+1kβkkk2NiΩkViWΩk+∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+μkkk2NiΩkDiWΩk+θkkk2NiΩkWiWΩk+vi+ζ8,i+1kβkkk2NiΩkViWΩk+∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+μkkk2NiΩkDiWΩk+θkkk2NiΩkWiWΩk=−ζ3,i+1k1−∑j=1pβkjViWΩj+μkjDiWΩj+θkjWiWΩjNiΩk+NiΩj−dk−vi+ζ2,i+1kβkkViWΩk2NiΩk+ζ4,i+1kμkkDiWΩk2NiΩk+ζ6,i+1kθkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pβkjkViWΩj+μkjkDiWΩj+θkjkWiWΩjNiΩk+NiΩj+βkkkViWΩk+μkkkDiWΩk+θkkkWiWΩk2NiΩk−ζ7,i+1kvi,(35)Δζ4,ik=−∂H∂DiWΩk=−A2−ζ1,i+1kγkk2NiΩkViMΩk−ζ3,i+1kμkk2NiΩkDiMΩk+ζ4,i+11−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩk+NiΩj−dk−vi−ζ5,i+1kνkk2NiΩkWiMΩk+ζ7,i+1kγkk2NiΩkViMΩk+μkkk2NiΩkDiMΩk+∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩk+NiΩj+νkkk2NiΩkWiMΩk+ζ8,i+1kγkk2NiΩkViMΩk+μkkk2NiΩkDiMΩk+∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩk+NiΩj+νkkk2NiΩkWiMΩk+vi=−A2−ζ4,i+1k1−∑j=1pγjkViMΩj+μjkDiMΩj+νjkWiMΩjNiΩj+NiΩk−dk−vi+ζ1,i+1kγkkViMΩk2NiΩk+ζ3,i+1kμkkDiMΩk2NiΩk+ζ5,i+1kνkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pγjkkViMΩj+μjkkDiMΩj+νjkkWiMΩjNiΩj+NiΩk+γkkkViMΩk+μkkkDiMΩk+νkkkWiMΩk2NiΩk−ζ8,i+1kvi,(36)Δζ5,ik=−∂H∂WiMΩk=−−ζ2,i+1kηkk2NiΩkViWΩk−ζ4,i+1kνkk2NiΩkDiWΩk+ζ5,i+11−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩj−dk−ζ6,i+1kσkk2NiΩkWiWΩk+ζ7,i+1kηkk2NiΩkViWΩk+νkkk2NiΩkDiWΩk+∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+σkkk2NiΩkWiWΩk+ζ8,i+1kηkk2NiΩkViWΩk+νkkk2NiΩkDiWΩk+∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+σkkk2NiΩkWiWΩk=−ζ5,i+1k1−∑j=1pηkjViWΩj+νkjDiWΩj+σkjWiWΩjNiΩk+NiΩj−dk+ζ2,i+1kηkkViWΩk2NiΩk+ζ4,i+1kνkkDiWΩk2NiΩk+ζ6,i+1kσkkWiWΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pηkjkViWΩj+νkjkDiWΩj+σkjkWiWΩjNiΩk+NiΩj+ηkkkViWΩk+νkkkDiWΩk+σkkkWiWΩk2NiΩk,(37)Δζ6,ik=−∂H∂WiWΩk=−−ζ1,i+1kδkk2NiΩkViMΩk−ζ3,i+1kθkk2NiΩkDiMΩk−ζ5,i+1kσkk2NiΩkWiMΩk+ζ6,i+11−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩk+NiΩj−dk+ζ7,i+1kδkk2NiΩkViMΩk+θkkk2NiΩkDiMΩk+σkkk2NiΩkWiMΩk+∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩk+NiΩj+ζ8,i+1kδkk2NiΩkViMΩk+θkkk2NiΩkDiMΩk+σkkk2NiΩkWiMΩk+∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩk+NiΩj=−ζ6,i+1k1−∑j=1pδjkViMΩj+θjkDiMΩj+σjkWiMΩjNiΩj+NiΩk−dk+ζ1,i+1kδkkViMΩk2NiΩk+ζ3,i+1kθkkDiMΩk2NiΩk+ζ5,i+1kσkkWiMΩk2NiΩk−ζ7,i+1k+ζ8,i+1k∑j=1pδjkkViMΩj+θjkkDiMΩj+σjkkWiMΩjNiΩj+NiΩk+δkkkViMΩk+θkkkDiMΩk+σkkkWiMΩk2NiΩk,(38)Δζ7,ik=−∂H∂MiMΩk=−ζ3,i+1kλk+ζ6,i+1kρkk+ζ7,i+1k1−λk+ρk,(39)Δζ8,ik=−∂H∂MiWΩk=−−A3+ζ4,i+1kλk+ζ5,i+1kωk+ζ8,i+1k1−λk+ωk,with transversality conditions:(40)ζ1,N=0,ζ2,N=A1,ζ3,N=0,ζ4,N=A2,ζ5,N=0,ζ6,N=0,ζ7,N=0ζ8,N=−A3.
To obtain the optimality conditions, we take the variation with respect to controlui∗,vi∗ and set it equal to zero:(41)∂H∂ui=τ1ui−ζ1,i+1kViMΩk−ζ2,i+1kViWΩk+ζ7,i+1kViMΩk+ζ8,i+1kViWΩk=0,∂H∂vi=τ2ui−ζ3,i+1kDiMΩk−ζ4,i+1kDiWΩk+ζ7,i+1kDiMΩk+ζ8,i+1kDiWΩk=0.
Then, we obtain the optimal control:(42)ui∗=ViMΩkζ1,i+1k−ζ7,i+1k+ViWΩkζ2,i+1k−ζ8,i+1kτ1;vi∗=DiMΩkζ3,i+1k−ζ7,i+1k+DiWΩkζ4,i+1k−ζ8,i+1kτ2.
By the bounds inUad, it is easy to obtain ui∗,vi∗ in the following form:(43)ui∗=minmaxViMΩkζ1,i+1k−ζ7,i+1k+ViWΩkζ2,i+1k−ζ8,i+1kτ1,umin,umax,vi∗=minmaxDiMΩkζ3,i+1k−ζ7,i+1k+DiWΩkζ4,i+1k−ζ8,i+1kτ2,umin,umax,for i=0,…,N−1.


## 4. Discussion

In this section, we provide numerical simulations to demonstrate our theoretical results in the case when the studied domain represents the assembly ofp regions (cities, towns, etc.). The code is written and compiled in MATLAB using the data cited in Table 2. The optimality systems are solved using an iterative method where, at instant i, the states ViM,ViW,…,MiW with an initial guess are obtained based on a progressive scheme in time, and their adjoint variables ζl,i+1k,l=1,2,…,8 are obtained based on a regressive scheme in time because of the transversality conditions. Afterward, we update the optimal control values (30) and (31) using the values of state and costate variables obtained in the previous steps. Finally, we execute the previous steps until a tolerance criterion is reached. In order to show the importance of our work and without loss of generality, we consider here that p=4. An area with 4 regions was considered: two regions Ω1 and Ω4 represent two urban areas and the regions Ω2 and Ω3 represent the two rural areas. The High Commissioner for Planning (HCP, public establishment) [14] gave various population statistics according to the regions (urban and rural). The values defined in Table 2 are inspired by data mentioned on page 2. Generally, the HCP gives the average sum of marriage, divorce, and widow rates between different population areas αkj, βkj, …, ωk. And, precisely, one of the main objects of our work is to estimate as close as possible the values of the parameters αkj, …, ωk defined in Table 2.Table 2
Parameters values associated with discrete-time systems (9)–(16) and with optimal control problem (18).
V0M=4500,V0W=4600,D0M=400D0W=600,W0M=200,W0W=300M0M=2500,M0W=2000Λ11=250,Λ12=250,d1=0.03Ω1:α1j11≤j≤4=0.07,0.05,0.06,0.05α1jj1≤j≤4=0.06,0.07,0.06,0.08β1j11≤j≤4=0.04,0.02,0.03,0.03β1jj1≤j≤4=0.03,0.04,0.04,0.05γ1j11≤j≤4=0.02,0.02,0.02,0.03γ1jj1≤j≤4=0.02,0.02,0.01,0.02δ1j11≤j≤4=0.01,0.02,0.01,0.01δ1jj1≤j≤4=0.01,0.01,0.01,0.01η1j11≤j≤4=0.03,0.02,0.02,0.03η1jj1≤j≤4=0.02,0.03,0.02,0.02μ1j11≤j≤4=0.02,0.02,0.01,0.02μ1jj1≤j≤4=0.01,0.02,0.01,0.01θ1j11≤j≤4=0.02,0.02,0.01,0.02θ1jj1≤j≤4=0.01,0.02,0.01,0.01ν1j11≤j≤4=0.02,0.01,0.01,0.02ν1jj1≤j≤4=0.01,0.02,0.01,0.01σ1j11≤j≤4=0.01,0.01,0.01,0.01σ1jj1≤j≤4=0.01,0.01,0.01,0.02λ1=0.003,0.001,0.001,0.004ρ1=0.01,0.01,0.01,0.01,ω1=0.03V0M=1050,V0W=1060,D0M=90D0W=100,W0M=40,W0W=45M0M=510,M0W=500Λ21=60,Λ22=65,d2=0.03Ω2:α2j21≤j≤4=0.06,0.07,0.04,0.06α2jj1≤j≤4=0.06,0.06,0.07,0.06β2j21≤j≤4=0.03,0.04,0.03,0.03β2jj1≤j≤4=0.03,0.03,0.04,0.04γ2j21≤j≤4=0.02,0.02,0.02,0.02γ2jj1≤j≤4=0.02,0.01,0.02,0.01δ2j21≤j≤4=0.01,0.02,0.01,0.01δ2jj1≤j≤4=0.01,0.01,0.02,0.01η2j21≤j≤4=0.02,0.02,0.03,0.02η2jj1≤j≤4=0.03,0.02,0.02,0.04μ2j21≤j≤4=0.01,0.02,0.02,0.01μ2jj1≤j≤4=0.01,0.02,0.01,0.02θ2j21≤j≤4=0.01,0.02,0.02,0.03θ2jj1≤j≤4=0.02,0.02,0.01,0.01ν2j21≤j≤4=0.01,0.02,0.02,0.01ν2jj1≤j≤4=0.01,0.01,0.01,0.02σ2j21≤j≤4=0.01,0.02,0.01,0.01σ2jj1≤j≤4=0.01,0.01,0.02,0.01λ2=0.002,0.001,0.001,0.004ρ2=0.01,0.01,0.01,0.01,ω2=0.02V0M=850,V0W=860,D0M=90D0W=110,W0M=40,W0W=50M0M=450,M0W=400Λ31=70,Λ32=75,d3=0.03Ω3:α3j11≤j≤4=0.05,0.06,0.05,0.05α3jj1≤j≤4=0.07,0.06,0.07,0.07β3j31≤j≤4=0.03,0.04,0.03,0.04β3jj1≤j≤4=0.03,0.02,0.04,0.03γ3j31≤j≤4=0.02,0.02,0.02,0.01γ3jj1≤j≤4=0.02,0.01,0.02,0.02δ3j31≤j≤4=0.01,0.01,0.01,0.02δ3jj1≤j≤4=0.01,0.02,0.01,0.01η3j31≤j≤4=0.04,0.03,0.03,0.02η3jj1≤j≤4=0.01,0.02,0.03,0.03μ3j31≤j≤4=0.02,0.01,0.03,0.02μ3jj1≤j≤4=0.01,0.02,0.01,0.01θ3j31≤j≤4=0.01,0.01,0.01,0.01θ3jj1≤j≤4=0.01,0.02,0.03,0.02ν3j31≤j≤4=0.01,0.01,0.01,0.02ν3jj1≤j≤4=0.01,0.02,0.01,0.01σ3j31≤j≤4=0.01,0.01,0.01,0.01σ3jj1≤j≤4=0.01,0.01,0.01,0.02λ3=0.01,0.02,0.01,0.01ρ3=0.01,0.01,0.01,0.01,ω3=0.03V0M=5500,V0W=5600,D0M=900D0W=820,W0M=400,W0W=300M0M=3500,M0W=3000Λ41=250,Λ42=260,d3=0.04Ω4:α4j41≤j≤4=0.05,0.04,0.06,0.07α4jj1≤j≤4=0.07,0.07,0.06,0.05β4j41≤j≤4=0.05,0.04,0.04,0.05β4jj1≤j≤4=0.03,0.02,0.03,0.01γ4j41≤j≤4=0.02,0.02,0.02,0.02γ4jj1≤j≤4=0.02,0.02,0.01,0.01δ4j41≤j≤4=0.01,0.01,0.02,0.01δ4jj1≤j≤4=0.01,0.02,0.01,0.01η4j41≤j≤4=0.03,0.03,0.02,0.04η4jj1≤j≤4=0.03,0.02,0.03,0.02μ4j41≤j≤4=0.02,0.01,0.02,0.01μ4jj1≤j≤4=0.01,0.02,0.01,0.01θ4j41≤j≤4=0.02,0.01,0.02,0.01θ4jj1≤j≤4=0.01,0.02,0.01,0.01ν4j41≤j≤4=0.02,0.01,0.02,0.01ν4jj1≤j≤4=0.01,0.02,0.01,0.01σ4j41≤j≤4=0.01,0.01,0.02,0.01σ4jj1≤j≤4=0.01,0.02,0.01,0.01λ4=0.001,0.001,0.003,0.003ρ4=0.01,0.01,0.01,0.01,ω4=0.03τ1=2.5×106,τ2=3×105A1=0.0005,A2=0.0004,A3=0.0002


### 4.1. Simulations without Controls

In this section, Figure1 depicts dynamics of the states ViM,ViW,…,MiW in regions Ω1, Ω2, Ω3, and Ω4, respectively, in the case when there is yet no control strategy, and we note that in all these figures presented here, simulations give us an idea about the evolution of marital status in each region whether rural or urban and the impact of flows of populations between regions on the dynamics of the marital status of each one.Figure 1
Time evolution of marital status of family dynamicsVMDW without controls.
(a)(b)(c)(d)The study of the evolution of the dynamics of the marital status was spread out over a period of 20 years. In Figure1, we note that the evolution of the population in regions Ω1 and Ω4 is almost the same and similarly for regions Ω2 and Ω3.In the urban areaΩ1, the number of single men is 4500 at time i=0, and the number of married men is about 2000 and 400 cases of divorce for men and the same numbers for women, with a birth rate of 250 per year and 0.03% mortality. The number of virgins declined a little and stabilized after 10 years to 3750. By cons, for the number of married men and women, there is a slight evolution that reaches after 12 years the number of 3200 married people and finally, for the number of divorces, it remains almost stable throughout the period and remains around 600 divorces.In ruralΩ2, the number of unmarried men at time i=0 is 1050, the number of married men is about 500 and 100 divorces for men and the same numbers for women, with a birth rate of 60 per year and 0.03% mortality. These numbers will be evolved, so the number of virgins has decreased to reach about 700 virgins. The number of married women will experience considerable growth to reach around 1400 for men and 1100 for women. The number of women divorces will witness a slight growth of about 300 divorces and stability for men of almost 100.


### 4.2. Simulations with Controls

Considering the critical level of control, we give optimal control sufficient to reduce the number of virgin and divorced individuals and to increase the married individuals inΩ1 and we also studied the effect of controls on other neighboring regions Ω2, Ω3, and Ω4.In the following, we can see that the optimal control function has a very desirable effect upon the population of virgin and divorced people which decreases while the married population increases in a consistent way during the length of the process. The time evolution of the respective populations with control is displayed in Figure2.Figure 2
Time evolution of marital status of family dynamicsVMDW with controls.
(a)(b)(c)(d)Figures3–5 allow us to compare changes in the number of virgins, divorced, and married individuals before and after the introduction of control.Figure 3
Time evolution of virgin and divorced populations with and without controls.
(a)(b)(c)(d)Figure 4
Time evolution of widowed and married populations with and without controls.
(a)(b)(c)(d)Figure 5
Time evolution of optimal controlsu∗ and v∗.
(a)(b)In a population of virgin men of 4500 individuals (Figure3(a)) by applying the control law, this number decreases rapidly to reach after 3 years the number of 1700 virgins and then goes back up to reach out at the end of the companion awareness 3000 virgin men. This result is similar to virgin women as shown in Figure 3(b).Figure3(c) also shows the effect of applying the control law by indicating that the number of divorced men decreases more rapidly at the beginning of treatment to reach 100 male divorces. Then, we notice that the number increases and we can justify it by the fact that the number of divorced people is proportional to the number of married people. Since the number of married individuals increases along the course of treatment, it has been assumed that divorced men remain in Ω1. Unlike divorced women who can return to their region of origin, therefore the number of divorced women that was 600 women decreases rapidly to reach 150 divorced women after 4 years and then increases to 280 divorced women at the end of the companion as shown in Figure 3(d).Figures4(c) and 4(d) show that the number of married men and women grows after 2 years and reaches a maximum value of about 2500 for married men and 2000 for married women; then it declines to 3800 for men and 4250 for women at the end of companion. This decline is due to the fact that the number of divorced men and women increases, but it is still quite high as far as the uncontrolled number of married men and women. In Figures 4(a) and 4(b), we show the evolution of widowed men and women; the number increases slightly concerning the number of widowed individuals in a dynamics of the marital system without control strategies.We also note in Figures1 and 2 and contrary to the evolution of the dynamics of marital status in Ω1 that the number of virgin and divorced people slightly increases in the regions Ω2, Ω3, and Ω4 and the number of married people decreases slightly.Finally, Figure5 displays the time evolution of the optimal control u∗,v∗.


## 4.1. Simulations without Controls

In this section, Figure1 depicts dynamics of the states ViM,ViW,…,MiW in regions Ω1, Ω2, Ω3, and Ω4, respectively, in the case when there is yet no control strategy, and we note that in all these figures presented here, simulations give us an idea about the evolution of marital status in each region whether rural or urban and the impact of flows of populations between regions on the dynamics of the marital status of each one.Figure 1
Time evolution of marital status of family dynamicsVMDW without controls.
(a)(b)(c)(d)The study of the evolution of the dynamics of the marital status was spread out over a period of 20 years. In Figure1, we note that the evolution of the population in regions Ω1 and Ω4 is almost the same and similarly for regions Ω2 and Ω3.In the urban areaΩ1, the number of single men is 4500 at time i=0, and the number of married men is about 2000 and 400 cases of divorce for men and the same numbers for women, with a birth rate of 250 per year and 0.03% mortality. The number of virgins declined a little and stabilized after 10 years to 3750. By cons, for the number of married men and women, there is a slight evolution that reaches after 12 years the number of 3200 married people and finally, for the number of divorces, it remains almost stable throughout the period and remains around 600 divorces.In ruralΩ2, the number of unmarried men at time i=0 is 1050, the number of married men is about 500 and 100 divorces for men and the same numbers for women, with a birth rate of 60 per year and 0.03% mortality. These numbers will be evolved, so the number of virgins has decreased to reach about 700 virgins. The number of married women will experience considerable growth to reach around 1400 for men and 1100 for women. The number of women divorces will witness a slight growth of about 300 divorces and stability for men of almost 100.


## 4.2. Simulations with Controls

Considering the critical level of control, we give optimal control sufficient to reduce the number of virgin and divorced individuals and to increase the married individuals inΩ1 and we also studied the effect of controls on other neighboring regions Ω2, Ω3, and Ω4.In the following, we can see that the optimal control function has a very desirable effect upon the population of virgin and divorced people which decreases while the married population increases in a consistent way during the length of the process. The time evolution of the respective populations with control is displayed in Figure2.Figure 2
Time evolution of marital status of family dynamicsVMDW with controls.
(a)(b)(c)(d)Figures3–5 allow us to compare changes in the number of virgins, divorced, and married individuals before and after the introduction of control.Figure 3
Time evolution of virgin and divorced populations with and without controls.
(a)(b)(c)(d)Figure 4
Time evolution of widowed and married populations with and without controls.
(a)(b)(c)(d)Figure 5
Time evolution of optimal controlsu∗ and v∗.
(a)(b)In a population of virgin men of 4500 individuals (Figure3(a)) by applying the control law, this number decreases rapidly to reach after 3 years the number of 1700 virgins and then goes back up to reach out at the end of the companion awareness 3000 virgin men. This result is similar to virgin women as shown in Figure 3(b).Figure3(c) also shows the effect of applying the control law by indicating that the number of divorced men decreases more rapidly at the beginning of treatment to reach 100 male divorces. Then, we notice that the number increases and we can justify it by the fact that the number of divorced people is proportional to the number of married people. Since the number of married individuals increases along the course of treatment, it has been assumed that divorced men remain in Ω1. Unlike divorced women who can return to their region of origin, therefore the number of divorced women that was 600 women decreases rapidly to reach 150 divorced women after 4 years and then increases to 280 divorced women at the end of the companion as shown in Figure 3(d).Figures4(c) and 4(d) show that the number of married men and women grows after 2 years and reaches a maximum value of about 2500 for married men and 2000 for married women; then it declines to 3800 for men and 4250 for women at the end of companion. This decline is due to the fact that the number of divorced men and women increases, but it is still quite high as far as the uncontrolled number of married men and women. In Figures 4(a) and 4(b), we show the evolution of widowed men and women; the number increases slightly concerning the number of widowed individuals in a dynamics of the marital system without control strategies.We also note in Figures1 and 2 and contrary to the evolution of the dynamics of marital status in Ω1 that the number of virgin and divorced people slightly increases in the regions Ω2, Ω3, and Ω4 and the number of married people decreases slightly.Finally, Figure5 displays the time evolution of the optimal control u∗,v∗.


## 5. Conclusions

In this paper, we consider a discrete-time marital status model; we treat the modelling and control of a system that described the case of a multiregion model. Optimal control is investigated to reduce the number of virgins and the divorced population and to increase the married ones. The first control introduced is supposing the benefits of an awareness campaign to educate virgin men and women about the benefits of marriage for individuals and the society, and the second control characterizes the legal procedures, administrative complications, and the heavy financial and social consequences of divorces. A discrete version of Pontryagin et al.’s maximum principle is done to analyze the optimal control problem and numerical simulation is given to illustrate the obtained results.



---

*Source: 9385261-2020-08-14.xml*