Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Modifying Models

Given a model that is already built and possibly optimized, the user might want to modify single constraint or variable bounds by means of correction or exploration of the feasible space.

In the following we show how single elements can be tweaked or rewritten. Let’s start with the simple model of the Getting Started section.

import pandas as pd
import xarray as xr

import linopy

m = linopy.Model()
time = pd.Index(range(10), name="time")

x = m.add_variables(
    lower=0,
    coords=[time],
    name="x",
)
y = m.add_variables(lower=0, coords=[time], name="y")

factor = pd.Series(time, index=time)

con1 = m.add_constraints(3 * x + 7 * y >= 10 * factor, name="con1")
con2 = m.add_constraints(5 * x + 2 * y >= 3 * factor, name="con2")

m.add_objective(x + 2 * y)
m.solve(solver_name="highs", output_flag=False)

m.solve(solver_name="highs", output_flag=False)
sol = m.solution.to_dataframe()
sol.plot(grid=True, ylabel="Optimal Value")

<Axes: xlabel='time', ylabel='Optimal Value'>

The figure above shows the optimal values of x(t) and y(t).

Varying lower and upper bounds

Now, let’s say we want to set the lower bound of x(t) to 1. This would translate to:

x.lower = 1

Note

The same could have been achieved by calling m.variables.x.lower = 1

Let’s solve it again!

m.solve(solver_name="highs", output_flag=False)
sol = m.solution.to_dataframe()
sol.plot(grid=True, ylabel="Optimal Value")

<Axes: xlabel='time', ylabel='Optimal Value'>

sol

	x	y
time
0	1.0	0.000000
1	1.0	1.000000
2	1.0	2.428571
3	1.0	3.857143
4	1.0	5.285714
5	1.0	6.714286
6	1.0	8.142857
7	1.0	9.571429
8	1.0	11.000000
9	1.0	12.428571

We see that the new lower bound of x is binding across all time steps.

Of course the implementation is flexible over the dimensions, so we can pass non-scalar values:

x.lower = xr.DataArray(range(10, 0, -1), coords=(time,))

m.solve(solver_name="highs", output_flag=False)
sol = m.solution.to_dataframe()
sol.plot(grid=True, ylabel="Optimal Value")

<Axes: xlabel='time', ylabel='Optimal Value'>

You can manipulate the upper bound of a variable in the same way.

Varying Constraints

A similar functionality is implemented for constraints. Here we can modify the left-hand-side, the sign and the right-hand-side.

Assume we want to relax the right-hand-side of the first constraint con1 to 8 * factor. This would translate to:

con1.rhs = 8 * factor

Note

The same could have been achieved by calling m.constraints.con1.rhs = 8 * factor

Let’s solve it again!

m.solve(solver_name="highs", output_flag=False)
sol = m.solution.to_dataframe()
sol.plot(grid=True, ylabel="Optimal Value")

<Axes: xlabel='time', ylabel='Optimal Value'>

In contrast to previous figure, we now see that the optimal value of y does not reach values above 10 in the end.

In the same way, we can modify the left-hand-side. Assume we want to weight y with a coefficient of 8 in the constraints, this gives

con1.lhs = 3 * x + 8 * y

Note: The same could have been achieved by calling

m.constraints['con1'].lhs = 3 * x + 8 * y

which leads to

m.solve(solver_name="highs", output_flag=False)
sol = m.solution.to_dataframe()
sol.plot(grid=True, ylabel="Optimal Value")

<Axes: xlabel='time', ylabel='Optimal Value'>

Varying the objective

Varying the objective happens in the same way as for the left-hand-side of the constraint as it is a linear expression too. Note, when passing an unstacked linear expression, i.e. an expression with more than the _term dimension, linopy will automatically stack it.

So assume, we would like to modify the weight of y in the objective function, this translates to:

m.objective = x + 3 * y

m.solve(solver_name="highs", output_flag=False)
sol = m.solution.to_dataframe()
sol.plot(grid=True, ylabel="Optimal Value")

<Axes: xlabel='time', ylabel='Optimal Value'>

As a consequence, y stays at zero for all time steps.

m.objective

Objective:
----------
LinearExpression: +1 x[0] + 3 y[0] + 1 x[1] ... +3 y[8] + 1 x[9] + 3 y[9]
Sense: min
Value: 139.0

Fixing Variables and Extracting MILP Duals

A common workflow in mixed-integer programming is to solve the MILP, then fix the integer/binary variables to their optimal values and re-solve as an LP to obtain dual values (shadow prices).

Let’s extend our model with a binary variable z that activates an additional capacity constraint on x.

z = m.add_variables(binary=True, coords=[time], name="z")

# x can only exceed 5 when z is active: x <= 5 + 100 * z
m.add_constraints(x <= 5 + 100 * z, name="capacity")

# Penalize activation of z in the objective
m.objective = x + 3 * y + 10 * z

m.solve(solver_name="highs", output_flag=False)

('ok', 'optimal')

Now fix the binary variable z to its optimal values and relax its integrality. This converts the model into an LP, which allows us to extract dual values.

m.variables.binaries.fix()
m.variables.binaries.relax()
m.solve(solver_name="highs", output_flag=False)

# Dual values are now available on the constraints
m.constraints["con1"].dual

<xarray.DataArray 'dual' (time: 10)> Size: 80B
array([-0.        , -0.        , -0.        ,  0.33333333,  0.33333333,
        0.375     ,  0.375     ,  0.375     ,  0.375     ,  0.375     ])
Coordinates:
  * time     (time) int64 80B 0 1 2 3 4 5 6 7 8 9

xarray.DataArray

'dual'

• time: 10

• -0.0 -0.0 -0.0 0.3333 0.3333 0.375 0.375 0.375 0.375 0.375

  array([-0.        , -0.        , -0.        ,  0.33333333,  0.33333333,
          0.375     ,  0.375     ,  0.375     ,  0.375     ,  0.375     ])

• Coordinates: (1)
  • time
    (time)
    int64
    0 1 2 3 4 5 6 7 8 9

    [10 values with dtype=int64]

Calling unfix() on all variables removes the fix constraints and unrelax() restores the integrality of z.

m.variables.unfix()
m.variables.unrelax()

# z is binary again
m.variables["z"].attrs["binary"]

True