The sensitivity toolbox provides a Pyomo interface to sIPOPT and k_aug to very quickly compute approximate solutions to nonlinear programs with a small perturbation in model parameters.
Using the Sensitivity Toolbox
We will start with a motivating example:
\[\begin{array}{ll}
\min_{x_1,x_2,x_3} & x_1^2 + x_2^2 + x_3^2 \\
\mathrm{s.t.} & 6 x_1 + 3 x_2 + 2 x_3 - p_1 = 0 \\
& p_2 x_1 + x_2 - x_3 - 1 = 0 \\
& x_1, x_2, x_3 \geq 0
\end{array}\]
Here \(x_1\), \(x_2\), and \(x_3\) are the decision variables while \(p_1\) and \(p_2\) are parameters. At first, let’s consider \(p_1 = 4.5\) and \(p_2 = 1.0\). Below is the model implemented in Pyomo.
# Import Pyomo and the sensitivity toolbox
>>> import pyomo.environ as pyo
>>> from pyomo.contrib.sensitivity_toolbox.sens import sensitivity_calculation
# Create a concrete model
>>> m = pyo.ConcreteModel()
# Define the variables with bounds and initial values
>>> m.x1 = pyo.Var(initialize = 0.15, within=pyo.NonNegativeReals)
>>> m.x2 = pyo.Var(initialize = 0.15, within=pyo.NonNegativeReals)
>>> m.x3 = pyo.Var(initialize = 0.0, within=pyo.NonNegativeReals)
# Define the parameters
>>> m.eta1 = pyo.Param(initialize=4.5,mutable=True)
>>> m.eta2 = pyo.Param(initialize=1.0,mutable=True)
# Define the constraints and objective
>>> m.const1 = pyo.Constraint(expr=6*m.x1+3*m.x2+2*m.x3-m.eta1 ==0)
>>> m.const2 = pyo.Constraint(expr=m.eta2*m.x1+m.x2-m.x3-1 ==0)
>>> m.cost = pyo.Objective(expr=m.x1**2+m.x2**2+m.x3**2)
The solution of this optimization problem is \(x_1^* = 0.5\), \(x_2^* = 0.5\), and \(x_3^* = 0.0\). But what if we change the parameter values to \(\hat{p}_1 = 4.0\) and \(\hat{p}_2 = 1.0\)? Is there a quick way to approximate the new solution \(\hat{x}_1^*\), \(\hat{x}_2^*\), and \(\hat{x}_3^*\)? Yes! This is the main functionality of sIPOPT and k_aug.
Next we define the perturbed parameter values \(\hat{p}_1\) and \(\hat{p}_2\):
>>> m.perturbed_eta1 = pyo.Param(initialize = 4.0)
>>> m.perturbed_eta2 = pyo.Param(initialize = 1.0)
And finally we call sIPOPT or k_aug:
>>> m_sipopt = sensitivity_calculation('sipopt', m, [m.eta1, m.eta2], [m.perturbed_eta1, m.perturbed_eta2], tee=False)
>>> m_kaug_dsdp = sensitivity_calculation('k_aug', m, [m.eta1, m.eta2], [m.perturbed_eta1, m.perturbed_eta2], tee=False)
The first argument specifies the method, either ‘sipopt’ or ‘k_aug’. The second argument is the Pyomo model. The third argument is a list of the original parameters. The fourth argument is a list of the perturbed parameters. It’s important that these two lists are the same length and in the same order.
First, we can inspect the initial point:
>>> print("eta1 = %0.3f" % m.eta1())
eta1 = 4.500
>>> print("eta2 = %0.3f" % m.eta2())
eta2 = 1.000
# Initial point (not feasible):
>>> print("Objective = %0.3f" % m.cost())
Objective = 0.045
>>> print("x1 = %0.3f" % m.x1())
x1 = 0.150
>>> print("x2 = %0.3f" % m.x2())
x2 = 0.150
>>> print("x3 = %0.3f" % m.x3())
x3 = 0.000
Next, we inspect the solution \(x_1^*\), \(x_2^*\), and \(x_3^*\):
# Solution with the original parameter values:
>>> print("Objective = %0.3f" % m_sipopt.cost())
Objective = 0.500
>>> print("x1 = %0.3f" % m_sipopt.x1())
x1 = 0.500
>>> print("x2 = %0.3f" % m_sipopt.x2())
x2 = 0.500
>>> print("x3 = %0.3f" % m_sipopt.x3())
x3 = 0.000
Note that k_aug does not save the solution with the original parameter values. Finally, we inspect the approximate solution \(\hat{x}_1^*\), \(\hat{x}_2^*\), and \(\hat{x}_3^*\):
# *sIPOPT*
# New parameter values:
>>> print("eta1 = %0.3f" %m_sipopt.perturbed_eta1())
eta1 = 4.000
>>> print("eta2 = %0.3f" % m_sipopt.perturbed_eta2())
eta2 = 1.000
# (Approximate) solution with the new parameter values:
>>> x1 = m_sipopt.sens_sol_state_1[m_sipopt.x1]
>>> x2 = m_sipopt.sens_sol_state_1[m_sipopt.x2]
>>> x3 = m_sipopt.sens_sol_state_1[m_sipopt.x3]
>>> print("Objective = %0.3f" % (x1**2 + x2**2 + x3**2))
Objective = 0.556
>>> print("x1 = %0.3f" % x1)
x1 = 0.333
>>> print("x2 = %0.3f" % x2)
x2 = 0.667
>>> print("x3 = %0.3f" % x3)
x3 = 0.000
# *k_aug*
# New parameter values:
>>> print("eta1 = %0.3f" %m_kaug_dsdp.perturbed_eta1())
eta1 = 4.000
>>> print("eta2 = %0.3f" % m_kaug_dsdp.perturbed_eta2())
eta2 = 1.000
# (Approximate) solution with the new parameter values:
>>> x1 = m_kaug_dsdp.x1()
>>> x2 = m_kaug_dsdp.x2()
>>> x3 = m_kaug_dsdp.x3()
>>> print("Objective = %0.3f" % (x1**2 + x2**2 + x3**2))
Objective = 0.556
>>> print("x1 = %0.3f" % x1)
x1 = 0.333
>>> print("x2 = %0.3f" % x2)
x2 = 0.667
>>> print("x3 = %0.3f" % x3)
x3 = 0.000