Frozen Models#

In rare cases, the GAMSPy model generation time dominates the solver solution time and GAMSPy itself becomes the bottleneck in an optimization application. If the model needs to be solved for multiple data sets (scenarios) the frozen model facility of GAMSPy can improve the overall performance. When one solve, GAMSPy generates the model instance from the model algebra with the current data, passes this to the solver, and after completion populates the variable and equation symbols in the container with the solution of the model instance. This model instance is created during the solve operation and discarded afterwards. This means that the next solve call will create a new model instance. The lifetime of a model instance of a frozen model goes beyond this and provides a controlled way of modifying a model instance and solving the resulting problem repeatedly in the most efficient way, by communicating whenever possible only the frozen model changes to the solver. Additionally, for continuous models, like LPs, some solvers can do a hot start from the previous basis solution without the need for refactorization. The terms “frozen model” and “model instance” are used interchangeably in the remainder of this text.

Frozen models are particularly useful when many scenarios of a model need to be solved where most of the data, especially the model rim, stays constant and only some (exogenous) parameters or bounds are subject to change. Typical examples are Monte-Carlo simulations with varying objective function coefficients.

The freeze method turns a regular model into a frozen model. This method actually generates/instantiates the model instance and requires in addition to model algebra all data of relevant model symbols. The freeze method gets a list of parameters and variable or equation attributes such as .lo, .up and .fx. Equation and variable attributes .l and .m in the list of modifiables are mainly used for starting non-linear models from different starting points.

The solve method takes the data of the modifiables from the container symbols and updates the frozen model. Variable and equation attributes are applied in a direct manner. Moreover, the data of the modifiable parameters is taken from the container symbol and applied to the model instance. Behind the scenes, the freeze method turned these modifiable parameters into fixed variables with the name of the parameter plus _var and mirrored them in the container. After a solve these variables hold the solution of the modifiable parameters. The marginal of these variables can provide useful sensitivity information about the parameter setting. This makes the _var suffix reserved: GAMS fixes every variable of the model instance with that suffix, therefore freeze rejects a model that contains one of its own. In addition, the container will contain the primal and dual solution with respect to the regular variables and equations of the model instance together with many model attributes including status, solve_status, objective_value, objective_estimation, solve_model_time, num_iterations, and num_nodes_used.

Note

Understanding the technical details that happen under the hood is not necessary for using frozen model effectively. Nevertheless, if you need to debug a frozen model, these details can help to interpret the representation of the model instance inside the solver. For non-linear models the fixed variables corresponding to the modifiables directly enter (non-linear) expressions in the model algebra. The automatic differentiation methods of GAMS and GAMSPy take these into account. For linear models, terms like a/sqr(b)*x, with the modifiable parameters a and b and variable x, become in a sense non-linear and need to be linearized. The Taylor expansion at the initial point allows such a linearization: \(f(x) = f(x.l) + \frac{d f(x.l)}{d x}(x-x.l)\). If f was linear with respect to the original variables, this is an exact representation and no error term is required. So our term a/sqr(b)*x will be represented (for better readability the a_var as been shortened to a, similar for b) as:

a.l/sqr(b.l)*x.l + a.l/sqr(b.l)*x + x.l/sqr(b.l)*a - 2*a.l*x.l/power(b.l,3)*b

or with an initial point of x.l=2 and fixed modifiers a.fx=6 and b.fx=3 we get

4/3 + 2/3*x + 2/9*a - 8/9*b

The following example demonstrates how to look at the solver representation of such a model instance:

import gamspy as gp

m = gp.Container()
x = gp.Variable(m, "x", records=2)
a = gp.Parameter(m, "a", records=6)
b = gp.Parameter(m, "b", records=3)
e = gp.Equation(m, "e", definition=a * x / gp.math.sqr(b) == 0)
mi = gp.Model(m, equations=[e], problem="LP", sense="FEASIBILITY")
mi.freeze([a, b])
mi.solve(solver="cplex", solver_options={"writelp": "mi.lp"})
mi.unfreeze()

The Cplex LP file mi.lp looks as follows:

\ENCODING=ISO-8859-1
\Problem name: gamsmodel

Minimize
 _obj: constobj
Subject To
 _e#0: 0.666666666666667 x + 0.222222222222222 a_var - 0.888888888888889 b_var
        = -1.33333333333333
Bounds
      x Free
      a_var = 6
      b_var = 3
      constobj = 0
End

In the solver’s presolve the fixed variables will be removed and there is minimal computational overhead with these newly introduced variables.

When all scenarios of a frozen model have been solved or the model needs to be changed beyond the modifiables, the unfreeze method needs to be called. This releases the resources of the model instance and turns the model back to a regular model that is generated and solved when the solve method is called. Any solution record that has not been read yet is read from the model instance before it is released, so the records of the symbols stay available afterwards.

Memory of a frozen model#

A frozen solve does not use the GAMS execution engine at all: the model, the solver and the data of the modifiables all live inside the Python process. The engine is only needed to generate the model instance, which happens in freeze. For a large model that matters, because the engine keeps the data it was given (and the model it generated) for as long as the Container is open, next to the model instance.

Passing hibernate=True to freeze puts the engine to sleep as soon as the model instance has been generated. The state of the engine is saved to a restart file before it is stopped, and the next statement that needs the engine restarts it from that state. Hibernation is therefore transparent: nothing has to be redeclared, the container stays usable, and the only thing that is given up in between is the memory the engine was holding:

model.freeze(modifiables=[ax], hibernate=True)   # engine stopped, its state saved

for scenario in scenarios:
    ax.setRecords(scenario)     # applied to the model instance, engine stays asleep
    model.solve()

How long the engine stays asleep depends on what the loop does. Setting the records of a modifiable does not synchronize with GAMS while the model is frozen, since a frozen solve reads them from the container, and the solve itself does not need the engine either. A statement that has to be executed by GAMS, on the other hand, such as ax[j] = 0 or any other assignment, wakes the engine up once, from the saved state, and it stays resident from then on. A loop that prepares its data with GAMS expressions has the engine back from its first iteration on; preparing that data in Python (setRecords) is what keeps it out of the picture for the whole scenario loop.

Saving a state and restarting from it is not free, so hibernating pays off only when the data the engine holds is large enough to be worth reclaiming. The same mechanism is available outside of freezing through Container.hibernate. It cannot be used while a Loop, For or While context manager is active, because the state of the engine cannot be saved before the loop is closed.

unfreeze reads the solution records that were never accessed before it releases the model instance, since they exist nowhere else. For a large model that is a lot of data, so unfreeze(load_records=False) discards them instead, and only the memory of the instance is left behind.

The following example shows how to use a single multiplier bmult to adjust the demand (equation demand) of the markets and solve the model with different value for this demand multiplier.

from gamspy import (
    Container,
    Set,
    Parameter,
    Variable,
    Equation,
    Sum,
    Model,
    Sense,
    ModelStatus,
)
import numpy as np

m = Container()

i = Set(m, name="i")
j = Set(m, name="j")

a = Parameter(
    m,
    name="a",
    domain=i,
    domain_forwarding=True,
    records=[["seattle", 350], ["san-diego", 600]],
)
b = Parameter(
    m,
    name="b",
    domain=j,
    domain_forwarding=True,
    records=[["new-york", 325], ["chicago", 300], ["topeka", 275]],
)
d = Parameter(
    m, name="d", domain=[i, j], records=np.array([[2.5, 1.7, 1.8], [2.5, 1.8, 1.4]])
)
c = Parameter(m, name="c", domain=[i, j])
c[i, j] = 90 * d[i, j] / 1000

x = Variable(m, name="x", domain=[i, j], type="Positive")

supply = Equation(m, name="supply", domain=i)
demand = Equation(m, name="demand", domain=j)
bmult = Parameter(m, name="bmult", records=1)

cost = Sum((i, j), c[i, j] * x[i, j])
supply[i] = Sum(j, x[i, j]) <= a[i]
demand[j] = Sum(i, x[i, j]) >= bmult * b[j]

transport = Model(
    m,
    name="transport",
    equations=[supply, demand],
    problem="LP",
    sense=Sense.MIN,
    objective=cost,
)

bmult_list = [0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3]

transport.freeze(modifiables=[bmult])

for b_value in bmult_list:
    bmult.setRecords(b_value)
    transport.solve(solver="conopt")
    print(
        f'obj:{transport.objective_value if transport.status == ModelStatus.OptimalGlobal else "infeasible"}'
    )

transport.unfreeze()

Note

  • Modifiable parameters cannot be used in .where conditions.

  • Variable and equation attributes used in equation algebra are evaluated once at model generation. Changes in the attributes will not percolate to the algebra. For example, the algebra x <= b * x.up will not change even if the modifiables include x.up. One needs a parameter bigM and algebra x <= b * bigM in order to modify this algebra in a frozen solve.

  • There are limitations to modifiable parameters in quadratic constraints in quadratic models (QCP, MIQCP, RMIQCP). For details, see Modifiables in Quadratic Models.

  • Even with few modifiable parameters, the frozen model may change significantly and the update of the frozen model may take longer than with a regular solve.

  • With a regular model the sparsity of the parameter often determines the sparsity of the constraint in a model. For example, Sum(j, a[j]*x[j]) will only have as many non-zeros as parameter a has. If you turn parameter a into a modifiable, this constraint will become totally dense, because we might update an element of a to a non-zero number even though the value is currently 0. One way around this is to limit the indexing of the Sum to relevant elements of j:

    jj = Set(m, "jj", domain=j)
jj[j] = a[j]
e = Equation(m, "e", definition=Sum(jj, a[jj]*x[jj] <= 1>))

  • Even though the solver can be switched in between solves, this is not recommended and can even lead to wrong results with FreezeOptions.update_type="accumulate" and solvers supporting communication of frozen model changes only.

Extra Values in the Modifiable’s Data#

The records of the modifiable parameters or attributes of a variable or equation might contain records that are not represented in the model. The FreezeOptions.no_match_limit allows to change the behavior in the presence of such extra records in the data.

In the following example, we define a standard LP model with coefficient matrix a, cost vector c, and modifiable right hand side vector b. The variable x is indexed over j, while the equations e are indexed over i. The actual model is defined over a subset of i and j, namely ii and jj:

import gamspy as gp

m = gp.Container()
i = gp.Set(m, name="i")
j = gp.Set(m, name="j")
a = gp.Parameter(m, name="a", domain=[i, j])
b = gp.Parameter(m, name="b", domain=i)
c = gp.Parameter(m, name="c", domain=j)

x = gp.Variable(m, name="x", domain=j, type="positive")
e = gp.Equation(m, name="e", domain=i)

ii = gp.Set(m, domain=i, name="ii", description="active i")
jj = gp.Set(m, domain=j, name="jj", description="active j")

e[ii] = gp.Sum(jj, a[ii, jj] * x[jj]) >= b[ii]

mymodel = gp.Model(
    m,
    name="mymodel",
    equations=[e],
    objective=gp.Sum(jj, c[jj] * x[jj]),
    sense="min",
    problem="lp",
)

i.setRecords(range(10))
j.setRecords(range(20))

ii.setRecords(range(5))
jj.setRecords(range(10))

a[ii, jj] = gp.math.uniform(0, 1)
b[ii] = gp.math.uniform(1, 10)
c[jj] = gp.math.uniform(1, 10)

mymodel.freeze(modifiables=[b])
b[i] = gp.math.uniform(1,10)  # assignment over the full set i, not over subset ii
mymodel.solve()

With this setup we can now pass more records to the modifiable parameter b than present in the frozen model, namely the b records 5, 6, 7, 8, and 9. If we have update values with no corresponding record in the model, this is suspicious and hence, the solve throws an exception by default:

gamspy.exceptions.GamspyException: Unmatched record limit exceeded while processing modifier b, for more info check no_match_limit option.

The number of allowed unmatched records can be controlled with no_match_limit. Here any number of 5 and larger will prevent the exception and the frozen model is passed to the solver:

mymodel.solve(freeze_options=gp.FreezeOptions(no_match_limit=5))

Missing Values in the Modifiable’s Data#

As extra values in the data cause concern, explicitly missing data for the parameters or variables present in the model also need some attention. GAMSPy and in particular the GAMS execution engine does not store default records. For parameters this means that zeros (0) are not explicitly stored. But also variable records with default attributes, e.g. positive variables with level, marginal, and lower bound 0.0, upper at +infinity, and scale at 1.0, are not explicitly stored. So GAMSPy does not distinguish between a missing or default value. If the frozen model contains a particular record of a modifiable symbol, but the corresponding symbol data in the container does not, GAMSPy treats this missing record to be at default value. There are several choices for how this default value should be applied when updating the frozen model. The different values for FreezeOptions.update_type control the behavior of default records in the modifiable’s data.

Basically, there are two options: Setting update_type="0" takes the values of the modifier’s data (default or not) and applies them to the records in the model. For parameters, this means, all records in the model are set to 0 (zero) and then the non-default values in the modifiable’s data are applied. For variables, this means the default bounds are restored and the records for the bounds stored in the variable symbols’s data are applied afterwards. Using update_type="base_case" or update_type="accumulate" will just ignore default values and only apply the explicit records found in the modifiable’s data. The option update_type="base_case" will apply the records against the model records present when the freeze was called; update_type="accumulate" will apply the records against the model records from the last solve. The difference will become clearer with the following example that builds on the example for no_match_limit. In order to make the updated b clearer we start the freeze with b=[1,2,3,4,5]:

b[ii] = gp.Ord(ii)
mymodel.freeze(modifiables=[b])
mymodel.solve()

Now, we clear b and set the last record to 100:

b[ii] = 0
b['4'] = 100

If we solve now with update_type="0", the solver sees the right hand side (which is our b parameter) of [0,0,0,0,100]. When update_type="base_case" or update_type="accumulate" is used, the solver sees the right hand side of [1,2,3,4,100]. Let’s continue the example with update_type="0":

mymodel.solve(freeze_options=gp.FreezeOptions(update_type="0"))

The records for b inside the model instance are now [0,0,0,0,100]. We could have also accomplished this by not having any default records and any update_type:

b[ii] = gp.SpecialValues.EPS
b['4'] = 100

Now assume, we again clear the updater symbol’s data and set the first record to 200

b[ii] = 0
b['0'] = 200

With update_type="0", the solver sees b=[200,0,0,0,0]; with update_type="base_case", the solver sees b=[200,2,3,4,5] and with update_type="accumulate", the solver sees b=[200,0,0,0,100] (assuming we solved the first scenario with update_type="0"). The default for this option is update_type="base_case".

The following example builds upon the model with the demand multiplier earlier in this chapter and shows how to update the upper bound of the transportation variable x to prevent any transportation on a particular connection:

transport.freeze(modifiables=[x.up])

for irec in i.toList():
    for jrec in j.toList():
        x.up[irec, jrec] = SpecialValues.EPS
        transport.solve()
        x.up[irec, jrec] = SpecialValues.POSINF # restore original bound
        print(
            f'obj:{transport.objective_value if transport.status == ModelStatus.OptimalGlobal else "infeasible"}'
        )

transport.unfreeze()

Modifiables in Quadratic Models#

Quadratic models have certain limitations when turned into frozen models. These limitations only affect solvers that extract the quadratic matrices from the non-linear expressions like Cplex, COPT, Gurobi, Mosek, and Xpress. Other solvers that work directly with the non-linear expressions and automatic differentiation like Conopt and Knitro do not have any limitations.

Variable and equation attributes, like variable bounds do not represent any limitation even for quadratic models. Only parameter modifiables can cause some problems. Since modifiable parameters are turned into fixed variables, any expression that is not linear or quadratic, like our example a/sqr(b)*x, are rejected with the message Detected n general nonlinear rows in model. Acceptance of terms that would be linear or constant in a regular model, but quadratic in a frozen model, like a*x or even a*b with modifiable parameters a and b and variable x depends heavily on the solver. Additive use of modifiable parameters, e.g. to change the right hand side of a constraint works without problems.

Due to limited update possibilities in the API of some solvers, quadratic models do not benefit from the minimal communication of just the model instance changes. The entire instance needs to be freshly set up inside the solver with every solve. Nevertheless, using a frozen model avoids regenerating the model instance from scratch for every solve, leading to improved efficiency.