Model

To create a model, you must instantiate the gamspy.Model class. This instance defines the scope of the problem, including which equations are included, the direction of optimization, and the problem type.

Essential Components

When initializing a model, you typically provide:

• Equations: A list of equations to define the constraints.

• Problem: The type of mathematical problem (e.g. LP, NLP, etc.).

• Sense: The direction of optimization (MIN, MAX, or FEASIBILITY).

• Objective: The expression or variable to be optimized.

Example: A Simple Knapsack Model

Here is the general syntax for defining a model, illustrated by a knapsack problem:

import gamspy as gp

# 1. Define Container and Symbols
m = gp.Container()
i = gp.Set(m, description="items")
p = gp.Parameter(m, description="profits", domain=i)
w = gp.Parameter(m, description="weights", domain=i)
c = gp.Parameter(m, description="capacity")
x = gp.Variable(m, domain=i, type=gp.VariableType.BINARY)

# 2. Define Equations
capacity_restriction = gp.Equation(m, definition=gp.Sum(i, w[i] * x[i]) <= c)

# 3. Instantiate the Model
knapsack = gp.Model(
    m,
    equations=m.getEquations(), # Automatically grabs all equations in the container
    problem=gp.Problem.MIP,     # Mixed Integer Program
    sense=gp.Sense.MAX,         # Maximize profit
    objective=gp.Sum(i, p[i] * x[i]),
)

Note

In case the name is not provided, it will be autogenerated.

Classification of Models

You must specify the correct model type (e.g., LP, NLP) before solving. GAMSPy verifies that your formulation matches the specified type and will issue an error if there is a mismatch (e.g., using nonlinear terms in a generic LP model).

Model Type	Model Type Description
LP	Linear Program
NLP	Nonlinear Program
QCP	Quadratically Constrained Program
DNLP	Discontinuous Nonlinear Program
MIP	Mixed Integer Program
RMIP	Relaxed Mixed Integer Program
MINLP	Mixed Integer Nonlinear Program
RMINLP	Relaxed Mixed Integer Nonlinear Program
MIQCP	Mixed Integer Quadratically Constrained Program
RMIQCP	Relaxed Mixed Integer Quadratically Constrained Program
MCP	Mixed Complementarity Problem
CNS	Constrained Nonlinear System
MPEC	Mathematical Programs with Equilibrium Constraints
RMPEC	Relaxed Mathematical Program with Equilibrium Constraints
EMP	Extended Mathematical Program
MPSGE	General Equilibrium

Tip

Some problems allow multiple valid types. For example, a model with binary variables can be solved as a MIP (exact integer solution) or an RMIP (Relaxed MIP, treating integers as continuous).

All model types are exposed with Problem enum, but the problem type can be specified as a string as well.

Also the direction types of the optimization (MIN, MAX, or FEASIBILITY) are exposed with Sense enum but it can be specified as a string similarly.

Solving a Model

Once defined, you trigger the optimization using the solve method.

Basic Solve

In most cases, calling model.solve() without arguments is sufficient. GAMSPy will use the default solver, problem type, sense and options for your problem type:

summary = model.solve()

The function returns a Pandas DataFrame summarizing the run, including the name of the solver, solver status, objective value, and solve time.

For example, the aformentioned knapsack problem would look like the following with the solve statement:

import gamspy as gp

m = gp.Container()
i = gp.Set(m, description="items")
p = gp.Parameter(m, description="profits", domain=i)
w = gp.Parameter(m, description="weights", domain=i)
c = gp.Parameter(m, description="capacity")
x = gp.Variable(m, "x", domain=i, type=gp.VariableType.BINARY)

capacity_restriction = gp.Equation(m, definition=gp.Sum(i, w[i] * x[i]) <= c)

knapsack = gp.Model(
    m,
    equations=m.getEquations(),
    problem=gp.Problem.MIP,
    sense=gp.Sense.MAX,
    objective=gp.Sum(i, p[i] * x[i]),
)
summary = knapsack.solve()  # This is the only line added.

Note

The data for the model was omitted on purpose to show that GAMSPy does not need data to define the model. See the full implementation with example data here: GAMS-dev/gamspy-examples

Specifying A Solver

In most cases, calling the solve function without any parameters is sufficient. In this scenario, the default solver depending on the problem type, default options will be used. However, users who require more control can specify the solver, general options, and solver-specific options.

All installed solvers on your system can be queried by running the following command:

gamspy list solvers

To see all available solvers that can be installed and used, run the following command.:

gamspy list solvers -a

Let’s say we want to use [HIGHS](https://highs.dev/) solver to solve the same knapsack problem. One can install highs with:

gamspy install solver highs

Then, change the solver as follow:

summary = knapsack.solve(solver="highs")  # This is the only line changed.

This solve statement will use highs solver instead of the default solver.

Redirecting Output

By default, solver output does not stream into the terminal to avoid cluttering your terminal. The output of the solver can be redirected to a file, to standard output or to any custom stream that supports write and flush operations by specifying the output parameter in the solve function.:

import gamspy as gp

m = gp.Container()
i = gp.Set(m, description="items")
p = gp.Parameter(m, description="profits", domain=i)
w = gp.Parameter(m, description="weights", domain=i)
c = gp.Parameter(m, description="capacity")
x = gp.Variable(m, "x", domain=i, type=gp.VariableType.BINARY)

capacity_restriction = gp.Equation(m, definition=gp.Sum(i, w[i] * x[i]) <= c)

knapsack = gp.Model(
    m,
    equations=m.getEquations(),
    problem=gp.Problem.MIP,
    sense=gp.Sense.MAX,
    objective=gp.Sum(i, p[i] * x[i]),
)

# redirect output to standard output (your console)
summary = knapsack.solve(output=sys.stdout)

# redirect output to a file
with open("my_out_file", "w") as file:
    knapsack.solve(output=file)

# redirect to custom stream
class MyStream:
    def write(self, data):
        logger.info(data.strip())

    def flush(self): ...

my_stream = MyStream()
knapsack.solve(output=my_stream)

Solver Options

You can fine-tune the solution process using Solve Options (generic GAMS options) and Solver Options (specific to the algorithm, like CPLEX or Gurobi settings).

Generic Solve Options

Solve options can be specified using the gamspy.Options() class. These control high-level behavior, such as iteration limits or time limits. For example:

import gamspy as gp

m = gp.Container()
... # Definition of your model
model = gp.Model(m, equations=m.getEquations(), problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=z)
model.solve(options=gp.Options(iteration_limit=2))

Commonly Used Options

• time_limit: Wall-clock time limit.

• iteration_limit: Maximum number of solver iterations.

• absolute_optimality_gap: Stop when the gap is below this threshold.

• threads: Number of threads/cores to use (if not defined, the solver decides on number of processors to use).

Here is the list of options and their descriptions:

Option	Description	Possible Values
cns	Default cns solver	Any solver installed in your system that can solve cns
dnlp	Default dnlp solver	Any solver installed in your system that can solve dnlp
emp	Default emp solver	Any solver installed in your system that can solve emp
lp	Default lp solver	Any solver installed in your system that can solve lp
mcp	Default mcp solver	Any solver installed in your system that can solve mcp
minlp	Default minlp solver	Any solver installed in your system that can solve minlp
mip	Default mip solver	Any solver installed in your system that can solve mip
miqcp	Default miqcp solver	Any solver installed in your system that can solve miqcp
mpec	Default mpec solver	Any solver installed in your system that can solve mpec
nlp	Default nlp solver	Any solver installed in your system that can solve nlp
qcp	Default qcp solver	Any solver installed in your system that can solve qcp
rminlp	Default rminlp solver	Any solver installed in your system that can solve rminlp
rmip	Default rmip solver	Any solver installed in your system that can solve rmip
rmiqcp	Default rmiqcp solver	Any solver installed in your system that can solve rmiqcp
rmpec	Default rmpec solver	Any solver installed in your system that can solve rmpec
license	Path to the license file	str
allow_suffix_in_equation	Allow variables with suffixes in model algebra	bool
allow_suffix_in_limited_variables	Allow domain limited variables with suffixes in model	bool
basis_detection_threshold	Basis detection threshold	float
compile_error_limit	Compile time error limit	int
domain_violation_limit	Domain violation limit solver default	int
generate_name_dict	Makes names of variables and equations that have been generated by the solve statement available to the solver	bool
enable_scaling	Determines whether to employ user-specified variable and equation scaling factors	bool
enable_prior	Instructs the solver to use the priority branching information. If and how priorities are used is solver-dependent.	bool
bypass_solver	If True, GAMSPy does not pass the generated model to the solver.	bool
hold_fixed_variables	Treat fixed variables as constants	bool
iteration_limit	Iteration limit of solver	int
keep_temporary_files	Controls keeping or deletion of process directory and scratch files	bool
listing_file	Listing file name	Name of the listing file
log_file	Log file name	Name of the log file
variable_listing_limit	Maximum number of columns listed in one variable block	int
equation_listing_limit	Maximum number of rows listed in one equation block	int
node_limit	Node limit in branch and bound tree	int
absolute_optimality_gap	Absolute Optimality criterion solver default	float
relative_optimality_gap	Relative Optimality criterion solver default	float
memory_tick_interval	Wait interval between memory monitor checks: ticks = milliseconds	float
monitor_process_tree_memory	Monitor the memory used by the GAMS process tree	bool
profile	Execution profiling	0: No profiling 1: Minimum profiling 2: Profiling depth for nested control structures
profile_file	Write profile information to this file	str
profile_tolerance	Minimum time a statement must use to appear in profile generated output	float
reference_file	Symbol reference file	str
time_limit	Wall-clock time limit for solver	float
savepoint	Save solver point in GDX file	0: No point GDX file is to be saved 1: A point GDX file from the last solve is to be saved 2: A point GDX file from every solve is to be saved 3: A point GDX file from the last solve is to be saved 4: A point GDX file from every solve is to be saved
seed	Random number seed	int
report_solution	Solution report print option	0: Remove solution listings following solves 1: Include solution listings following solves 2: Suppress all solution information
show_os_memory		0: Show memory reported by internal accounting 1: Show resident set size reported by operating system 2: Show virtual set size reported by operating system
solve_link_type	Solve link option	“disk”: The model instance is saved to the scratch directory, “memory”: The model instance is passed to the solver in-memory
merge_strategy	Multiple solve management	“replace” | “merge” | “clear”
step_summary	Summary of computing resources used by job steps	bool
suppress_compiler_listing	Compiler listing option	bool
report_solver_status	Solver Status file reporting option	bool
threads	Number of threads to be used by a solver	int
write_listing_file	Controls listing file creation	bool
zero_rounding_threshold	The results of certain operations will be set to zero if abs(result) LE ZeroRes	float
report_underflow	Report underflow as a warning when abs(results) LE ZeroRes and result set to zero	bool

To check all available options, see gamspy.Options().

Note

Solve options can also be created with a dictionary of GAMS options. For example:

import gamspy as gp

m = gp.Container()
... # Definition of your model
model = gp.Model(m, equations=m.getEquations(), problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=z)
model.solve(options=gp.Options.fromGams({"reslim": 5, "lp": "gurobi"}))

Solver-Specific Options

In addition to solve options, user can specify solver options as a dictionary. These are passed as a dictionary and vary by solver (e.g., Gurobi, CONOPT).:

import gamspy as gp

m = gp.Container()
... # Definition of your model
model = gp.Model(m, equations=m.getEquations(), problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=z)
model.solve(solver="conopt", solver_options={"rtmaxv": "1.e12"})

For all possible solver options, please check the corresponding solver manual.

Solver options can also be created with gamspy.Container.writeSolverOptions():

import gamspy as gp

m = gp.Container()
...
...
...
your model definition goes here
...
...
...
m.writeSolverOptions(solver="nlpec", solver_options={"testTol": 1e-006})
model.solve(solver="reshop", options=Options(mcp="nlpec"), solver_options={"subsolveropt": 1})

In this example, gamspy.Container.writeSolverOptions() would create solver options for the nlpec solver.

Execution Backends

GAMSPy allows you to solve models locally or on cloud platforms.

Solving Locally

By default, models are solved locally (on your machine).

Solving with GAMS Engine

Synchronous Solve

In order to send your model to be solved with GAMS Engine, you need to define the GAMS Engine configuration. This can be done by importing EngineClient and creating an instance. The user can then pass this instance to the solve method and specify the backend as engine.

import gamspy as gp

m = gp.Container()
z = gp.Variable(m, description="objective variable")
e1 = gp.Equation(m)
e1[...] = ... # definition of the equation
e2 = gp.Equation(m)
e2[...] = ... # definition of the equation

model = gp.Model(m, equations=[e1, e2], problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=z)

client = gp.EngineClient(
    host=os.environ["ENGINE_URL"],
    username=os.environ["ENGINE_USER"],
    password=os.environ["ENGINE_PASSWORD"],
    namespace=os.environ["ENGINE_NAMESPACE"],
)
model.solve(solver="conopt", backend="engine", client=client)

Asynchronous Solve

If you just want to send your jobs to GAMS Engine without blocking until the results are received, the is_blocking parameter can be set to False in EngineClient.

Tokens of the submitted jobs are stored in client.tokens

client = gp.EngineClient(
    host=os.environ["ENGINE_URL"],
    username=os.environ["ENGINE_USER"],
    password=os.environ["ENGINE_PASSWORD"],
    namespace=os.environ["ENGINE_NAMESPACE"],
    is_blocking=False,
)

for _ in range(3):
    ... # changes in your model
    model.solve(backend="engine", client=client)

print(client.tokens) # This prints all tokens for the submitted jobs

The results of the non-blocking jobs can be retrieved later. For example if want to retrieve the results of the last submitted job, we can do that following:

token = client.tokens[-1]
client.job.get_results(token, working_directory="out_dir")

The results would be downloaded to the given working directory. The downloaded GDX file will have the same name with gdxOutputPath. Then, if one wants to read the results, they can simply create a new Container and read the results from the downloaded GDX file:

gdx_out_path = os.path.join("out_dir", os.path.basename(m.gdxOutputPath()))
solution_container = gp.Container(load_from=gdx_out_path)

Solving with NEOS Server

Synchronous Solve

In order to send your model to be solved to NEOS Server, you need to create a NeosClient. This can be done by importing NeosClient and creating an instance. The user can then pass this instance to the solve method and specify the backend as neos.

import gamspy as gp

m = gp.Container()
z = gp.Variable(m, description="objective variable")
e1 = gp.Equation(m)
e1[...] = ... # definition of the equation
e2 = gp.Equation(m)
e2[...] = ... # definition of the equation

model = gp.Model(m, equations=[e1, e2], problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=z)

client = gp.NeosClient(
    email=os.environ["NEOS_EMAIL"],
    username=os.environ["NEOS_USER"],
    password=os.environ["NEOS_PASSWORD"],
)
model.solve(backend="neos", client=client)

Providing your username and password is optional for the NEOS Server backend, but it is recommended as it allows you to review your models on the NEOS web client. The environment variables can be set in a .env file or with export statements on the command line. Example of running your model on NEOS Server without authentication:

NEOS_EMAIL=<your_email> python <your_script>

If one wants to investigate the results later on NEOS Server web client, they can provide the username and password in the same way:

NEOS_EMAIL=<your_email> NEOS_USER=<your_username> NEOS_PASSWORD=<your_password> python <your_script>

Alternatively, the output of NEOS can be redirected to a file by specifying the output stream:

model.solve(backend="neos", client=client, output=sys.stdout)

Note

NEOS Server backend does not support loadpoint option and external equations at the moment.

Asynchronous Solve

If you just want to send your jobs to NEOS server without blocking until the results are received, is_blocking parameter can be set to False in NeosClient.

All submitted jobs are stored in client.jobs in case you want to reach to the job numbers and job passwords you already sent to the server.

client = gp.NeosClient(
    email=os.environ["NEOS_EMAIL"],
    username=os.environ["NEOS_USER"],
    password=os.environ["NEOS_PASSWORD"],
    is_blocking=False,
)

for _ in range(3):
    ... # changes in your model
    model.solve(backend="neos", client=client)

print(client.jobs) # This prints all job numbers and jon passwords as a list of tuples

The results of the non-blocking jobs can be retrieved later. For example if want to retrieve the results of the last submitted job, we can do that following:

job_number, job_password = client.jobs[-1]
client.get_final_results(job_number, job_password)
client.download_output(job_number, job_password, working_directory="my_out_directory")

The results would be downloaded to the given working directory. The downloaded gdx file will always have the name “output.gdx”. Then, if one wants to read the results, they can simply create a new Container and read the results from the downloaded gdx file:

solution_container = gp.Container(load_from="my_out_directory/output.gdx")

The terms of use for NEOS can be found here: Terms of use.

Model Attributes

Models have attributes that store a variety of information, including

• information about the results of solving a model, the results of a solve, and the model’s solution,

• information about certain features to be used by GAMSPy or the solver,

• information passed to GAMSPy or the solver specifying various settings that are also available as option.

Model Attribute	Description
num_domain_violations	Number of domain violations
algorithm_time	Solver-dependent timing information
total_solve_time	Elapsed time it took to execute a solve statement in total
total_solver_time	Elapsed time taken by the solver only
num_iterations	Number of iterations used
marginals	Indicator for marginals present
max_infeasibility	Maximum of infeasibilities
mean_infeasibility	Mean of infeasibilities
status	Integer number that indicates the model status
num_nodes_used	Number of nodes used by the MIP solver
solve_number	Number of the last solve
num_dependencies	Number of dependencies in a CNS model
num_discrete_variables	Number of discrete variables
num_equations	Number of equations
num_infeasibilities	Number of infeasibilities
num_nonlinear_insts	Number of nonlinear instructions
num_nonlinear_zeros	Number of nonlinear nonzeros
num_nonoptimalities	Number of nonoptimalities
num_nonzeros	Number of nonzero entries in the model coefficient matrix
num_mcp_redefinitions	Number of MCP redefinitions
num_variables	Number of variables
num_bound_projections	Number of bound projections during model generation
objective_estimation	Estimate of the best possible solution for a mixed-integer model
objective_value	Objective function value
used_model_type	Integer number that indicates the used model type
model_generation_time	Time GAMS took to generate the model in wall-clock seconds
solve_model_time	Time the solver used to solve the model in seconds
sum_infeasibilities	Sum of infeasibilities
solve_status	Indicates the solver termination condition
solver_version	Solver version

Converting A Model To A Scalar Format

gamspy.Model.convert() transforms a GAMSPy model instance into a scalar model where all confidential information has been removed or into formats used by other modeling and solution systems. It is designed to achieve the following goals:

• Permit users to convert a confidential model into GAMS scalar format so that any idenifiable structure is removed. It can then be passed on to others for investigation without confidentiality being lost.

• A way of sharing a model instance created by GAMSPy for use with other modeling systems or solvers.

For example, you can convert your GAMSPy model into a scalar GAMS model as follows:

import gamspy as gp

m = gp.Container()
... # Definition of your model
model = gp.Model(m, equations=m.getEquations(), problem="LP", sense="MAX", objective=z)
model.convert(path="path_to_the_directory", file_format=gp.FileFormat.GAMS, options=gp.ConvertOptions(Width=50))

The path parameter specifies the directory where the converted model will be saved. The file_format parameter specifies the format of the converted model. The options parameter specifies the options for the conversion. See gamspy.ConvertOptions() for all available conversion options.

Exporting Model To LaTeX

GAMSPy models can be exported to a .tex file in LaTeX format by using the toLatex function of the model. The generated .tex file can be automatically compiled into a PDF file by using xelatex

import gamspy as gp

m = gp.Container()
... # Definition of your model
model = gp.Model(m, equations=m.getEquations(), problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=z)
model.toLatex(path=<latex_path>, generate_pdf=True)

Note

CMU Serif font must be installed in your machine to see Greek letters in your pdf document.

Note

To generate a PDF file from a .tex file, you must install xelatex on your system and add it to your PATH. The toLatex function uses the names of the GAMSPy symbols. If names are not supplied, GAMSPy invents (ugly) names which would show up in the LaTeX source. So for this feature to be useful the GAMSPy set, parameter, variable, and equations should be specified with a name.

GAMSPy also supports renaming symbols for LaTeX output in the toLatex function. This can be done by providing a dictionary mapping the GAMSPy symbol names to LaTeX names. For example:

import gamspy as gp

m = gp.Container()
... # Definition of your model
variance = gp.Variable(m)
model = gp.Model(m, equations=m.getEquations(), problem=gp.Problem.LP, sense=gp.Sense.MAX, objective=variance)

rename = {"variance": r"\ensuremath{\sigma}"}

model.toLatex(path=<latex_path>, rename=rename)

Latex representation of the individual equation definitions can be retrieved with equation.latexRepr. For example:

import gamspy as gp

m = gp.Container()
i = gp.Set(m)
v = gp.Variable(m, domain=i)
e = gp.Equation(m, domain=i)
e[i] = v[i] >= 5
print(e.latexRepr())

This would result in:

$v_{i} \geq 5\hfill \forall i$

Matches for MCP Models

Mixed Complementarity Problem (MCP) models can be defined as pair-wise complementarities between variables and equations. The Model class accepts these pair-wise complementarities via the matches argument in its constructor.

p = gp.Variable(m, type=gp.VariableType.POSITIVE, domain=c)
y = gp.Variable(m, type=gp.VariableType.POSITIVE, domain=s)
i = gp.Variable(m, type=gp.VariableType.POSITIVE, domain=h)

mkt = gp.Equation(m, domain=c)
profit = gp.Equation(m, domain=s)
income = gp.Equation(m, domain=h)

mkt[c] = gp.Sum(s, a[c, s] * y[s]) + gp.Sum(h, e[c, h]) >=
         gp.Sum(h.where[esub[h] != 1],
                (i[h] / gp.Sum(cc, alpha[cc, h] * p[cc] ** (1 - esub[h])))
              * alpha[c, h]
              * (1 / p[c]) ** esub[h],
            ) + gp.Sum(h.where[esub[h] == 1], i[h] * alpha[c, h] / p[c])
profit[s] = -gp.Sum(c, a[c, s] * p[c]) >= 0
income[h] = i[h] >= gp.Sum(c, p[c] * e[c, h])

hansen = gp.Model(
    m,
    problem=gp.Problem.MCP,
    matches={mkt: p, profit: y, income: i},
)

You do not need to include equations already provided in matches in the equations argument.

In addition to this explicit equation, variable matching, some alternative matching constructs with more flexibility are also supported.

Equation sequence syntax:

model = gp.Model(m, problem=gp.Problem.MCP, matches={(e1, e2, e3) : v})

This syntax requires that each equation in the sequence be conformant with v (i.e. each equation has the same domain with the variable) and that the set of tuples defining each equation be disjoint. For each column of v, at most one of e1 or e2 or e3 will have a matching row. This is useful when a variable contains columns of different kinds, e.g. prices for both tradable and non-tradable commodities whose equilibrium conditions are best expressed in different equations.

Variable sequence syntax:

model = gp.Model(m, problem=gp.Problem.MCP, matches={e : (v1, v2, v3)})

This construct requires that each variable in the variable sequence be conformant with e. This points to the exclusive-or relationship among the non-fixed variables involved in a match. For each row of e, at most one of the matching columns in v1 or v2 or v3 is allowed to be non-fixed. The fixed columns in the match are ignored by the solver, and the row is paired with the one non-fixed column. If all the columns are fixed, this effectively drops the row from the model, just as would happen with a fixed column in the simple match e.v. If no columns exist to match a row of e, this is an error. This construct is useful when a system has too many degrees of freedom if all the variables in question are left endogenous: by fixing some variables (i.e. making them exogenous) we arrive at a square system.

More background information on MCP models can be found here. An example MCP model can be found in the model library: HANSMCP.

Limiting Domain for Variables

It is possible to limit the domain of variables used in a model in the Model constructor. This allows to restrict the generation of blocks of variables in a single place instead of using, e.g., where statements at every place where this variable block is used in equations.

The following examples are based on the transportation model. To limit the transportation network in that model to certain links (e.g. because some are blocked because of some reason) one could introduce a subset of the possible links and use that with where conditions (equation supply) or indexing a subset (equation demand) in the equations like this:

import gamspy as gp

# Define symbols here
...

# Initialize whole network as free
free_links = gp.Set(
    m, domain=[i,j], description="usable links in the network",
    records=i.toList() + j.toList()
)
cost[...] = z == gp.Sum((i,j), (c[i,j] * x[i,j]).where[free_links[i,j]])
supply[i] = gp.Sum(j, x[i,j].where[free_links[i,j]]) <= a[i]
demand[j] = gp.Sum(free_links[i,j], x[i,j]) >= b[j]

# Block a particular link
free_links['san-diego','topeka'] = False

transport = gp.Model(
    m, equations=m.getEquations(), problem="LP",
    sense="MIN", objective=z
)
transport.solve()

Instead of adding the where condition or index subset to each appearance of x in the model, one could simply add a domain restriction for that variable to the model statement directly by specifying a variable and the set that limits its domain. Using this approach, the previous example looks like the following:

import gamspy as gp

# Define symbols here
...

# Initialize whole network as free
free_links = gp.Set(
    m, domain=[i,j], description="usable links in the network",
    records=i.toList() + j.toList()
)
cost[...] = z == gp.Sum((i,j), c[i,j] * x[i,j])
supply[i] = gp.Sum(j, x[i,j]) <= a[i]
demand[j] = gp.Sum(i, x[i,j]) >= b[j]

# Block a particular link
free_links['san-diego','topeka'] = False

transport = gp.Model(
    m, limited_variables=[x[free_links]], equations=m.getEquations(),
    problem="LP", sense="MIN", objective=z
)
transport.solve()