Two-Strain Tuberculosis Model#

Mathematical Formulation#

Problem Statement#

Find the optimal treatment controls \(u_1(t)\) and \(u_2(t)\) that minimize the total infected population and control costs:

\[J = \int_0^{t_f} \left[ L_2 + I_2 + \frac{1}{2}B_1 u_1^2 + \frac{1}{2}B_2 u_2^2 \right] dt\]

Subject to the epidemiological dynamics:

\[\frac{dS}{dt} = \Lambda - \beta_1 S \frac{I_1}{N} - \beta^* S \frac{I_2}{N} - \mu S\]
\[\frac{dT}{dt} = u_1 r_1 L_1 - \mu T + (1-(1-u_2)(p+q)) r_2 I_1 - \beta_2 T \frac{I_1}{N} - \beta^* T \frac{I_2}{N}\]
\[\frac{dL_1}{dt} = \beta_1 S \frac{I_1}{N} - (\mu + k_1) L_1 - u_1 r_1 L_1 + (1-u_2) p r_2 I_1 + \beta_2 T \frac{I_1}{N} - \beta^* L_1 \frac{I_2}{N}\]
\[\frac{dL_2}{dt} = (1-u_2) q r_2 I_1 - (\mu + k_2) L_2 + \beta^* (S + L_1 + T) \frac{I_2}{N}\]
\[\frac{dI_1}{dt} = k_1 L_1 - (\mu + d_1) I_1 - r_2 I_1\]
\[\frac{dI_2}{dt} = k_2 L_2 - (\mu + d_2) I_2\]

Boundary Conditions#

  • Initial conditions: \(S(0) = 19000\), \(T(0) = 250\), \(L_1(0) = 9000\), \(I_1(0) = 1000\), \(L_2(0) = 500\), \(I_2(0) = 250\)

  • Final conditions: All states free

  • Control bounds: \(0.05 \leq u_1, u_2 \leq 0.95\)

Physical Parameters#

  • \(\Lambda = 429\) per year

  • \(\beta_1 = 13\), \(\beta_2 = 13\), \(\beta^* = 0.029\) per year

  • \(\mu = 0.0143\) per year

  • \(d_1 = 0\), \(d_2 = 0\) per year

  • \(k_1 = 0.5\), \(k_2 = 1\) per year

  • \(r_1 = 2\), \(r_2 = 1\) per year

  • \(p = 0.4\), \(q = 0.1\)

  • \(N = 30000\)

  • \(B_1 = 50\), \(B_2 = 500\)

State Variables#

  • \(S(t)\)

  • \(T(t)\)

  • \(L_1(t)\)

  • \(I_1(t)\)

  • \(L_2(t)\)

  • \(I_2(t)\)

Control Variables#

  • \(u_1(t)\)

  • \(u_2(t)\)

Notes#

This model represents the optimal control of a two-strain tuberculosis epidemic with drug-sensitive and drug-resistant strains, where treatment decisions must balance infection reduction with control costs.

Running This Example#

cd examples/two-strain_tuberculosis_model
python two-strain_tuberculosis_model.py

Code Implementation#

examples/two-strain_tuberculosis_model/two-strain_tuberculosis_model.py#
  1import numpy as np
  2
  3import maptor as mtor
  4
  5
  6# Problem parameters from Table 10.30
  7beta1 = 13
  8beta2 = 13
  9mu = 0.0143
 10d1 = 0
 11d2 = 0
 12k1 = 0.5
 13k2 = 1
 14r1 = 2
 15r2 = 1
 16p = 0.4
 17q = 0.1
 18N = 30000
 19beta_star = 0.029
 20B1 = 50
 21B2 = 500
 22Lambda = mu * N  # 429
 23
 24# Problem setup
 25problem = mtor.Problem("Two-Strain Tuberculosis Model")
 26phase = problem.set_phase(1)
 27
 28# Variables
 29t = phase.time(initial=0.0, final=5.0)
 30S = phase.state("S", initial=19000.0)  # 76*N/120
 31T = phase.state("T", initial=250.0)  # N/120
 32L1 = phase.state("L1", initial=9000.0)  # 36*N/120
 33I1 = phase.state("I1", initial=1000.0)  # 4*N/120
 34L2 = phase.state("L2", initial=500.0)  # 2*N/120
 35I2 = phase.state("I2", initial=250.0)  # N/120
 36
 37u1 = phase.control("u1", boundary=(0.05, 0.95))
 38u2 = phase.control("u2", boundary=(0.05, 0.95))
 39
 40# Dynamics (equations 10.1009-10.1014)
 41phase.dynamics(
 42    {
 43        S: Lambda - beta1 * S * (I1 / N) - beta_star * S * (I2 / N) - mu * S,
 44        T: (
 45            u1 * r1 * L1
 46            - mu * T
 47            + (1 - (1 - u2) * (p + q)) * r2 * I1
 48            - beta2 * T * (I1 / N)
 49            - beta_star * T * (I2 / N)
 50        ),
 51        L1: (
 52            beta1 * S * (I1 / N)
 53            - (mu + k1) * L1
 54            - u1 * r1 * L1
 55            + (1 - u2) * p * r2 * I1
 56            + beta2 * T * (I1 / N)
 57            - beta_star * L1 * (I2 / N)
 58        ),
 59        L2: ((1 - u2) * q * r2 * I1 - (mu + k2) * L2 + beta_star * (S + L1 + T) * (I2 / N)),
 60        I1: k1 * L1 - (mu + d1) * I1 - r2 * I1,
 61        I2: k2 * L2 - (mu + d2) * I2,
 62    }
 63)
 64
 65# Objective: minimize L2 + I2 + 0.5*B1*u1^2 + 0.5*B2*u2^2
 66integrand = L2 + I2 + 0.5 * B1 * u1**2 + 0.5 * B2 * u2**2
 67integral_var = phase.add_integral(integrand)
 68problem.minimize(integral_var)
 69
 70# Mesh and guess
 71phase.mesh([6, 6, 6], [-1.0, -1 / 3, 1 / 3, 1.0])
 72
 73states_guess = []
 74controls_guess = []
 75for N_interval in [6, 6, 6]:
 76    tau = np.linspace(-1, 1, N_interval + 1)
 77    t_norm = (tau + 1) / 2
 78
 79    # Linear interpolation for states
 80    S_vals = 19000.0 - 1000.0 * t_norm
 81    T_vals = 250.0 + 50.0 * t_norm
 82    L1_vals = 9000.0 - 500.0 * t_norm
 83    I1_vals = 1000.0 - 200.0 * t_norm
 84    L2_vals = 500.0 - 100.0 * t_norm
 85    I2_vals = 250.0 - 50.0 * t_norm
 86
 87    states_guess.append(np.vstack([S_vals, T_vals, L1_vals, I1_vals, L2_vals, I2_vals]))
 88
 89    # Control guess (mid-range values)
 90    u1_vals = np.full(N_interval, 0.5)
 91    u2_vals = np.full(N_interval, 0.5)
 92    controls_guess.append(np.vstack([u1_vals, u2_vals]))
 93
 94phase.guess(states=states_guess, controls=controls_guess, integrals=5000.0)
 95
 96# Solve
 97solution = mtor.solve_adaptive(
 98    problem,
 99    error_tolerance=1e-6,
100    max_iterations=20,
101    min_polynomial_degree=3,
102    max_polynomial_degree=12,
103    nlp_options={
104        "ipopt.print_level": 0,
105        "ipopt.max_iter": 3000,
106        "ipopt.tol": 1e-8,
107        "ipopt.constr_viol_tol": 1e-7,
108    },
109)
110
111# Results
112if solution.status["success"]:
113    print(f"Objective: {solution.status['objective']:.5f}")
114    print(
115        f"Reference: 5152.07310 (Error: {abs(solution.status['objective'] - 5152.07310) / 5152.07310 * 100:.3f}%)"
116    )
117
118    # Final state values
119    S_final = solution[(1, "S")][-1]
120    T_final = solution[(1, "T")][-1]
121    L1_final = solution[(1, "L1")][-1]
122    I1_final = solution[(1, "I1")][-1]
123    L2_final = solution[(1, "L2")][-1]
124    I2_final = solution[(1, "I2")][-1]
125
126    print("Result:")
127    print(f"S: {S_final:.1f}")
128    print(f"T: {T_final:.1f}")
129    print(f"L1: {L1_final:.1f}")
130    print(f"I1: {I1_final:.1f}")
131    print(f"L2: {L2_final:.1f}")
132    print(f"I2: {I2_final:.1f}")
133
134    solution.plot()
135else:
136    print(f"Failed: {solution.status['message']}")