<jupyter_start><jupyter_text>ProblemWe want to maximize the total number of chairs and tables by using available wood and workers. |           | Wood | Worker ||-----------|------|--------|| Chair     | 1    | 2      || Table     | 2    | 1      || Available | 6    | 6      | The table show the resources required (on average) for an ideal case.  Ideal case - simple linear programming Assume the resource consumption is exact.<jupyter_code># Coefficients of the cost function
cost = -np.array([1., 1.]) # convert into a minimize problem

# Coefficients of the constraints
cons1 = np.array([1., 2.])
cons2 = np.array([2., 1.])

# Cost function
fun = lambda x: np.sum(cost * x)

# Constraints - ineq: greater than 0
cons = ({'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons1 * x)},
        {'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons2 * x)})

# Bounds
bnds = ((0, None),(0, None))

# Initial guess
guess = np.array([1., 1.])
# Optimize
res = opt.minimize(fun, guess, method = 'SLSQP', bounds = bnds, constraints = cons, tol = 1e-8)
res
res.x
print("objective:", np.sum(cost * res.x))
print("constraint 1:", round(np.sum(cons1 * res.x), 6) <= 6)
print("constraint 2:", round(np.sum(cons2 * res.x), 6) <= 6)<jupyter_output>objective: -4.0
constraint 1: True
constraint 2: True<jupyter_text>Deal with uncertainty What if the resource consumption is uncertaion? i.e., the coeffiecients deviate form its nominal (average) value.   |           | Wood    | Worker  ||-----------|---------|---------|| Chair     | 1 ± 0.8 | 2 ± 1   || Table     | 2 ± 1.5 | 1 ± 0.5 || Available | 6       | 6       |  Optmistic case We will use the least resource to produce more units.|           | Wood          | Worker        ||-----------|---------------|---------------|| Chair     | 1 - 0.8 = 0.2 | 2 - 1 = 1     || Table     | 2 - 1.5 = 0.5 | 1 - 0.5 = 0.5 || Available | 6             | 6             |<jupyter_code># Coefficients of the cost function
cost = -np.array([1., 1.]) # convert into a minimize problem

# Coefficients of the constraints
cons1 = np.array([0.2, 0.5])
cons2 = np.array([1.0, 0.5])

# Cost function
fun = lambda x: np.sum(cost * x)

# Constraints - ineq: greater than 0
cons = ({'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons1 * x)},
        {'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons2 * x)})

# Bounds
bnds = ((0, None),(0, None))

# Initial guess
guess = np.array([1., 1.])
# Optimize
res = opt.minimize(fun, guess, method = 'SLSQP', bounds = bnds, constraints = cons, tol = 1e-8)
res
res.x
print("objective:", np.sum(cost * res.x))
print("constraint 1:", round(np.sum(cons1 * res.x), 6) <= 6)
print("constraint 2:", round(np.sum(cons2 * res.x), 6) <= 6)<jupyter_output>objective: -12.0
constraint 1: True
constraint 2: True<jupyter_text>Pessimistic case We will use the most resource to produce less units.|           | Wood          | Worker        ||-----------|---------------|---------------|| Chair     | 1 + 0.8 = 1.8 | 2 + 1 = 3     || Table     | 2 + 1.5 = 3.5 | 1 + 0.5 = 1.5 || Available | 6             | 6             |<jupyter_code># Coefficients of the cost function
cost = -np.array([1., 1.]) # convert into a minimize problem

# Coefficients of the constraints
cons1 = np.array([1.8, 3.5])
cons2 = np.array([3.0, 1.5])

# Cost function
fun = lambda x: np.sum(cost * x)

# Constraints - ineq: greater than 0
cons = ({'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons1 * x)},
        {'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons2 * x)})

# Bounds
bnds = ((0, None),(0, None))

# Initial guess
guess = np.array([1., 1.])
# Optimize
res = opt.minimize(fun, guess, method = 'SLSQP', bounds = bnds, constraints = cons, tol = 1e-8)
res
res.x
print("objective:", np.sum(cost * res.x))
print("constraint 1:", round(np.sum(cons1 * res.x), 6) <= 6)
print("constraint 2:", round(np.sum(cons2 * res.x), 6) <= 6)<jupyter_output>objective: -2.46153846221
constraint 1: True
constraint 2: True<jupyter_text>Robust optimization The optmization model can represent as this form:\begin{array}{r l}    \max        & z = \sum_{j=1}^{n} c_j x_j & \\    \text{s.t.} & \sum_{j=1}^{n} a_{ij} x_{j} \leq b_i & \forall\, i=1,\dotsc,m \\                & x_j \geq 0 & \forall\, j=1,\dotsc,n \\\end{array} , where $a_{ij} \in [\bar a_{ij} -\hat a_{ij},\bar a_{ij} +\hat a_{ij}]$ We can define the **uncerrtainty set** as:$\mathcal U_i = \{ a_{ij}\in \mathbb {R} : a_{ij}= \bar a_{ij} + \hat a_{ij}d_{ij}, -1 \leq d_{ij} \leq 1, \sum_{j \in J_i} d_{ij} \leq \Gamma_i  \}$,where $\Gamma_i$ can be at most equal to the number of the uncertain parameters in the $i^{th}$ contraint; $J_i$ is the set of indexes of the constraints. \begin{array}{r l}\max & z = \sum_{j=1}^{n} c_j x_j & \\\text{s.t.} & \sum_{j=1}^{n} \bar a_{ij} x_{j} + \max_{S_i \subseteq J_i, |S_i| = \Gamma_i} { \sum_{j \in S_i} \hat a_{ij} y_j } \leq b_i & \forall\, i=1,\dotsc,m \\& x_j \leq y_j & \forall\, j=1,\dotsc,n\\& x_j \geq 0,y_j \geq 0  & \forall\, j=1,\dotsc,n\end{array} If $\Gamma_i = 0$, then it becomes the nominal problem; if m > $\Gamma_i$ > 0, we are looking for the worst-case. Use strong duality to take care of the max in the constraint, and substitute the the objective of the dual problem into the constraint.\begin{array}{r l}\max & z = \sum_{j=1}^{n} c_j x_j \\s.t.& \sum_{j=1}^{n} \bar a_{ij} x_{j} + z_i \Gamma_i + \sum_{j \in J_i} p_{ij} \leq b_i & \forall i \\& z_i + p_{ij} \geq \hat a_{ij} y_j & \forall i, j\in J_i\\& x_j \leq y_j & \forall j=1,...,n\\& p_{ij} \geq 0 & \forall i, j \in J_i\\& z_i \geq 0 & \forall i\\& x_j \geq 0,y_j \geq 0  & \forall j\end{array} Therefore, the original problem can be reformulated to:\begin{array}{r l}\max & z =x_1 + x_2 \\s.t.& 2x_1 + x_2 + z_1 \Gamma_1 + p_{11} + p_{12} \leq 6 \\& x_1 + 2x_2 + z_2 \Gamma_2 + p_{21} + p_{22} \leq 6 \\& z_1 + p_{11} \geq y_1\\& z_1 + p_{12} \geq 0.5 y_2\\& z_2 + p_{21} \geq 0.8 y_1\\& z_2 + p_{22} \geq 1.5 y_2\\& x_1 \leq y_1 \\& x_2 \leq y_2 \\& p_{ij} \geq 0 ,\forall i,j\\& z_1,z_2 \geq 0\\& y_1, y_2 \geq 0\\& x_1,x_2 \geq 0\end{array}<jupyter_code>def robust(g):
    # Coefficients of the cost function
    cost = -np.array([1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) # convert into a minimization problem

    # Coefficients of the constraints
    # [x1, x2, y1, y2, z1, z2, p11, p12, p21, p22]
    cons1 = np.array([ 2.0, 1.0, 0.0, 0.0,g[0], 0.0, 1.0, 1.0, 0.0, 0.0])
    cons2 = np.array([ 1.0, 2.0, 0.0, 0.0, 0.0,g[1], 0.0, 0.0, 1.0, 1.0])
    cons3 = np.array([ 0.0, 0.0,-1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0])
    cons4 = np.array([ 0.0, 0.0, 0.0,-0.5, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    cons5 = np.array([ 0.0, 0.0,-0.8, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0])
    cons6 = np.array([ 0.0, 0.0, 0.0,-1.5, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0])
    cons7 = np.array([-1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    cons8 = np.array([ 0.0,-1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

    # Cost function
    fun = lambda x: np.sum(cost * x)

    # Constraints - ineq: greater than 0
    cons = ({'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons1 * x)},
            {'type': 'ineq', 'fun': lambda x: 6 - np.sum(cons2 * x)},
            {'type': 'ineq', 'fun': lambda x: np.sum(cons3 * x)},
            {'type': 'ineq', 'fun': lambda x: np.sum(cons4 * x)},
            {'type': 'ineq', 'fun': lambda x: np.sum(cons5 * x)},
            {'type': 'ineq', 'fun': lambda x: np.sum(cons6 * x)},
            {'type': 'ineq', 'fun': lambda x: np.sum(cons7 * x)},
            {'type': 'ineq', 'fun': lambda x: np.sum(cons8 * x)}
           )

    # Bounds
    bnds = ((0, None),) * 10

    # Initial guess
    guess = np.ones(10)

    # Optimize
    res = opt.minimize(fun, guess, method = 'SLSQP', bounds = bnds, constraints = cons, tol = 1e-8)
    
    return res.x[0:2]<jupyter_output><empty_output><jupyter_text>If we asumme the same level of uncertainty, i.e., gamma_1 = gamma_2<jupyter_code>robust([0, 0])
x = np.linspace(0, 2, 21)
decision = [robust([i, i]) for i in x]
decision
bar1 = []
bar2 = []

for i in decision:
    bar1.append(i[0])
    bar2.append(i[1])

plt.figure(figsize=(10, 5))
plt.bar(x, bar1, width = 0.05, align='center')
plt.bar(x, bar2, width = 0.05, align='center', bottom = bar1)
plt.xlabel("Gamma_1 = Gamma_2")
plt.ylabel("Number of units")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2141309581126625, 'fraction_numerical': 0.0601829561868079, 'mean_word_length': 2.5799224472210254, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '383729', 'license': 'MIT', 'n_tokens_mistral': 3657, 'n_tokens_neox': 3366, 'n_words': 1012, 'path': 'maximize_units.ipynb', 'repo_name': 'jerrychens/robust_optimization', 'size': '8308'}