1import pyomo.environ as pe
2from pyomo.contrib import appsi
3from pyomo.common.timing import HierarchicalTimer
4
5
6def main(plot=True, n_points=200):
7    import numpy as np
8
9    # create a Pyomo model
10    m = pe.ConcreteModel()
11    m.x = pe.Var()
12    m.y = pe.Var()
13    m.p = pe.Param(initialize=1, mutable=True)
14
15    m.obj = pe.Objective(expr=m.x**2 + m.y**2)
16    m.c1 = pe.Constraint(expr=m.y >= (m.x + 1)**2)
17    m.c2 = pe.Constraint(expr=m.y >= (m.x - m.p)**2)
18
19    opt = appsi.solvers.Cplex()  # create an APPSI solver interface
20    opt.config.load_solution = False  # modify the config options
21    opt.update_config.check_for_new_or_removed_vars = False  # change how automatic updates are handled
22    opt.update_config.update_vars = False
23
24    # write a for loop to vary the value of parameter p from 1 to 10
25    p_values = [float(i) for i in np.linspace(1, 10, n_points)]
26    obj_values = list()
27    x_values = list()
28    timer = HierarchicalTimer()  # create a timer for some basic profiling
29    timer.start('p loop')
30    for p_val in p_values:
31        m.p.value = p_val
32        res = opt.solve(m, timer=timer)
33        assert res.termination_condition == appsi.base.TerminationCondition.optimal
34        obj_values.append(res.best_feasible_objective)
35        opt.load_vars([m.x])
36        x_values.append(m.x.value)
37    timer.stop('p loop')
38    print(timer)
39
40    if plot:
41        import matplotlib.pyplot as plt
42        # plot the results
43        fig, ax1 = plt.subplots()
44        ax1.set_xlabel('p')
45        ax1.set_ylabel('objective')
46        ax1.plot(p_values, obj_values, ':k', label='objective')
47
48        ax2 = ax1.twinx()
49        ax2.set_ylabel('x')
50        ax2.plot(p_values, x_values, '-b', label='x')
51
52        fig.legend()
53        plt.show()
54
55
56if __name__ == '__main__':
57    main()
58