An optimal control model of the spread of the COVID-19 pandemic in Iraq: Deterministic and chance-constrained model

3.1 Deterministic model

Optimal control for optimization is defined by: Min . J = ∑ t = 0 T { 1 + y ( t ) } 2(1)

Subject to the state equations of cases, inpatients, recovered cases and death: Δ C ( t ) = { 1 + y ( t ) } β ( t ) C ( t )(2)

Δ I ( t ) = { 1 + y ( t ) } β ( t ) C ( t ) - { μ ( t ) + N ( t ) } I ( t )(3) Δ R ( t ) = μ ( t ) I ( t )(4) Δ D ( t ) = N ( t ) I ( t )(5)

where,

y(t): The control variable.

C(t): Percent of confirmed cases.

I(t): Percent of inpatients.

R(t): Percent of recovered cases.

D(t): Percent of death cases.

β(t): Percent of the new cases to the inpatient cases.

μ(t): Percent of the recover cases to the inpatient cases.

N ( t ) : Percent of the death cases to the inpatient cases.

ΔC(t) = C(t)-C(t-1)

Equation (2) signifies the total number of cases that increased based on new cases daily. The total number of inpatients increases by the addition of new cases and decreases by the number of recovery and deaths each day (Equation 3). Equations 5) represent the total number of recovered patients and deaths, respectively. The control variable represents the percentage of isolation or each weather factor.

By using the PMP, determining the solution of the optimal control model can be achieved. The Lagrangian function is expressed as follows [17]: L = ∑ t = 0 T - { 1 + y ( t ) } 2 + ∑ t = 0 T λ c ( t ) [ { 1 + y ( t ) } β ( t ) C ( t ) - C ( t ) + C ( t - 1 ) ] + ∑ t = 0 T λ i ( t ) [ { 1 + y ( t ) } β ( t ) C ( t ) - { μ ( t ) + N ( t ) } I ( t ) - I ( t ) + I ( t - 1 ) ] + ∑ t = 0 T λ r ( t ) [ μ ( t ) I ( t ) - R ( t ) + R ( t - 1 ) ] + ∑ t = 0 T λ d ( t ) [ N ( t ) I ( t ) - D ( t ) + D ( t - 1 ) ](6)

A Hamiltonian function is expressed as: H ( t ) = - { 1 + y ( t ) } 2 + λ c ( t ) [ { 1 + y ( t ) } β ( t ) C ( t ) ] + λ i ( t ) [ { 1 + y ( t ) } β ( t ) C ( t ) - { μ ( t ) + N ( t ) } I ( t ) ] + λ r ( t ) [ μ ( t ) I ( t ) ] + λ d ( t ) [ N ( t ) I ( t ) ](7)

Substituting Equation (7) into Equation (6), gives: L = ∑ t = 0 T [ H ( t ) - λ c ( t ) { C ( t ) - C ( t - 1 ) } - λ i ( t ) { I ( t ) - I ( t - 1 ) } - λ r ( t ) { R ( t ) - R ( t - 1 ) } - λ d ( t ) { D ( t ) - D ( t - 1 ) } ](8)

Deriving Equation (8) concerning C (t) , I (t) , R (t) , D (t), separately, gives:

{ 1 + y ( t ) } β ( t ) { λ c ( t ) + λ i ( t ) } + λ c ( t ) - λ c ( t - 1 ) = 0 ,(9)

{ μ ( t ) + N ( t ) } λ i ( t ) - μ ( t ) λ r ( t ) - N ( t ) λ i ( t ) - λ i ( t ) + λ i ( t - 1 ) = 0 ,(10) λ r ( t ) - λ r ( t - 1 ) = 0 ,(11) λ d ( t ) - λ d ( t - 1 ) = 0 ,(12)

By rearranging Equations (9–12), we can determine the adjoint equations as: Δ λ c ( t ) = - { 1 + y ( t ) } β ( t ) { λ c ( t ) + λ i ( t ) }(13)

Δ λ i ( t ) = { μ ( t ) + N ( t ) } λ i ( t ) - μ ( t ) λ r ( t ) - N ( t ) λ d ( t )(14) Δ λ r ( t ) = 0 ,(15) Δ λ d ( t ) = 0 .(16)

From Equations 14), we get:

λ c ( t ) = 1 1 + { 1 + y ( t ) } β ( t ) [ λ c ( t - 1 ) + { 1 + y ( t ) } β ( t ) λ i ( t ) ] .(17)

λ i ( t ) = 1 1 - μ ( t ) - N ( t ) { λ i ( t - 1 ) - μ ( t ) λ r ( t ) - N ( t ) λ d ( t ) } .(18)

Rearranging Equation (13), yields: { 1 + y ( t ) } β ( t ) = - Δ λ c ( t ) { λ c ( t ) + λ i ( t ) }(19)

By substituting Equation (19) into Equation (2), it yields: Δ C ( t ) = - Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } C ( t )(20)

From Equation (20), we can determine the total number of cases over time as follows: C ( t ) = 1 1 - [ Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } ] C ( t - 1 )(21)

By substituting Equation (19) into Equation (3), it yields: Δ I ( t ) = - Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } C ( t ) - { μ ( t ) + N ( t ) } I ( t )(22)

From Equation (22), we then get: I ( t ) = 1 1 + μ ( t ) + N ( t ) [ I ( t - 1 ) - Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } C ( t ) ](23)

Then, substituting Equation (21) into Equation (23), it yields: I ( t ) = 1 1 + μ ( t ) + N ( t ) [ I ( t - 1 ) - { Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } } { 1 1 - [ Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } ] } C ( t - 1 ) ](24)

Equation (24) represents the total number of inpatients over time. We can determine the total number of recovered cases over time by substituting Equation (24) into Equation (4) as follows:

R ( t ) = R ( t - 1 ) + μ ( t ) 1 + μ ( t ) + N ( t ) [ I ( t - 1 ) - { Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } } { 1 1 - [ Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } ] } C ( t - 1 ) ](25)

Substituting Equation (24) into Equation (5) yields the total number of death cases over time as follows:

D ( t ) = D ( t - 1 ) + N ( t ) 1 + μ ( t ) + N ( t ) [ I ( t - 1 ) - { Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } } { 1 1 - [ Δ λ c ( t ) { λ c ( t ) + λ i ( t ) } ] } C ( t - 1 ) ](26)

Thus, the value of the control variable is as follows: y ( t ) = { - 1 ; no new cases ( optimal ) ( - 1 , 0 ) ; decrease cases 0 ; no effect > 0 ; increase cases(27)

From Equation (27), y (t) = -1 means: Δ λ c ( t ) = Δ C ( t ) = 0(28) Δ I ( t ) = - { μ ( t ) + N ( t ) } I ( t )(29) λ i ( T ) = 0(30)

Initially the value of the control variable is zero, which means the model represents the actual cases registered in Iraq with the real percentage of new cases, recovered, and deaths as follows: β 1 = C 1 - C 0 C 1 ; μ 1 = R 1 I 1 ; N 1 = D 1 I 1(31)

Next, by introducing the effect of wind (for example): changing the value of the control variable and cases depending on the wind degree (W) and transition matrix (φ) (see Appendix) we get: δ 1 = C 1 - C 0 + ABS ( W 1 - W 0 ) * φ(32) new C 1 = old C 0 + δ 1 ; new C 2 = new C 1 + δ 2(33) y 1 = δ 1 β 1 C 1 - 1(34)

Where ABS = absolute value.

The elements of the transition matrix can be calculated as follows (for example, humidity): φ ij = NC ¯ H i - NC ¯ H j H i - H j(35)

Where NC ¯ H i the average of new cases, according to the humidity degree H _i (Table 3 shows the average of new cases).

Finally, the next tasks include incorporating the value of the control variable into an optimal control model to obtain the solution. The solution of the optimal control model relies on Equations (15–18, 21, 24–26). The solution is found by using the goal seek function in Microsoft Excel with λ(T) = 0. Hence, achieving the condition λ(T) = 0 by changing the value of λ(0).

The number of new cases reported in Iraq depends on the number of samples tested. Therefore, the actual cases may be is greater than those recorded cases. In this model, the number of new cases is uncertain given by the equation representing chance-constrained: P r { Δ C ( t ) ⩾ NC } ⩾ α(36)

Where α takes values between zero and one (1) and NC is a random variable (new cases).

In determining the solution, chance-constrained must be converted to deterministic constrained. Here, the value of α is equal to zero or one (1) representing an extremely risky or extremely conservative attitude, respectively. The minimum acceptable to achieve the constraint is α, while (1 - α) is the maximum acceptable risks [7].

If NC is a random variable that adheres to a normal distribution with mean NC ¯ and standard deviation σ _NC , then the deterministic constraint that is equivalent to chance-constraint (36) as follows [4]: P r { Δ C ( t ) - NC ¯ σ NC ⩾ NC - NC ¯ σ NC } ⩾ α ∅ { Δ C ( t ) - NC ¯ σ NC } ⩾ α Δ C ( t ) ⩾ NC ¯ + σ NC ∅ - 1 ( α )(37)

Where ∅ is the cumulative distribution function (CDF) of the standard normal distribution.

The number of daily cases can be found from Eq. (37) with the initial value of cases C (0). To determine the effect of the weather factors, we apply the deterministic model (Equations 1–5) with the daily cases of chance-constrained (Equation 37).