@@ 1,186 @@
+# Xpress Demo
+
+This is a demo for the FICO Xpress module for linear programming optimization.
+Everything shown here is a 'normal' Python conversion of the
+[Jupyter notebook](https://gist.github.com/sburtner/3c058b06fc56bf5293c9851dd645939a)
+from the 2021 Spatial Data Science Symposium session
+"Spatial Optimization for Planning and Decision-Making".
+All credit and ownership belongs to Susan Burtner, Dr. Jing Xu, Jiwon Baik,
+Seonga Cho, Vanessa Echeverri Figueroa, Evgeny Noi, B. Amelia Pludow, and
+Enbo Zhou.
+
+## Relevant Notebook Quotes
+
+> Here, we are going to work on the [well-known] Knapsack problem, apply it to
+> forest treatment, and explore some extensions.
+
+> The context that we are interested in is to select land management units for
+> monitoring and evaluation in Stanislaus National Forest, in order to maximize
+> the total benefit from treatment while not exceeding the capacity of
+> resources and labor. The study area is part of Stanislaus National Forest
+> which is located northern California.
+
+> Mathematical programming usually has the following form:
+>
+> Maximize/Minimize f(X_1_, X_2_, ... x_n_)
+>
+> Subject to g_1_(X_1_, X_2_, ... x_n_) <= b_1_; g_2_(X_1_, X_2_, ... x_n_) <= b_2_; ... g_m_(X_1_, X_2_, ... x_n_) <= b_m_
+>
+> Where X_i_, i = 1, 2, ... n, are decision variables,
+> f is the objective function measuring soluton quality,
+> and we optimize f to acquire the best solution.
+> g_i_, i = 1, 2, ... m, are mathmatical [functions] to help formalize
+> constraints.
+> b_i_, i = 1, 2, ... m are constants to define constraints.
+
+> [L]et us assume that that the US Forest Service is planning to treat 100
+> acres (threshold capacity) of the Stanislaus National Forest to decrease
+> wildfire risk.
+> The selected units should be those that maximize the benefit or resilience
+> to wildfire while the threshold capacity is not surpassed.
+
+> Most of the real-world problems are based on multiple objectives.
+> ...
+> For example, (1) protecting property and (2) protecting people can be two
+> major objectives (Church and Murray, 2018) in fire protection services.
+> And most real-world problem contains a complicated relationship between
+> multiple objectives like the example.
+
+> The result of the multi-objective optimization problem is usually visualized
+> as a Pareto frontier.
+> The graph below shows that there is a trade-off relationship between
+> (1) covered property and (2) covered people.
+
+*Adjustments noted with brackets. Omissions noted with elipses.*
+
+
+## Results
+
+```
+$ make run
+(. .venv/bin/activate; python -m main)
+
+:: This is the data
+(16, 8)
+ CSO_depart CSOterr_ID CNSM_NoMGT Acres LMU_ID res_d_STND res_d_SPM geometry
+0 0.833334 164 172765.0 5.337479 35564 0.921725 14.163780 MULTIPOLYGON (((757755.000 4243875.000, 757725...
+1 0.565610 164 1757300.0 27.354568 36492 1.304450 20.044963 MULTIPOLYGON (((756915.000 4243875.000, 756975...
+2 0.585455 164 3106670.0 48.926867 36493 1.220160 18.749712 MULTIPOLYGON (((756735.000 4243545.000, 756765...
+3 0.590000 164 1042850.0 15.790036 36494 0.687985 10.571991 MULTIPOLYGON (((756465.000 4243695.000, 756435...
+4 0.652973 164 435944.0 8.228611 36964 1.695560 26.054994 MULTIPOLYGON (((757425.000 4243665.000, 757425...
+
+
+:: This is some licensing stuff, can be ignored
+Using the Community license in this session. If you have a full Xpress license, pass the full path to your license file to xpress.init(). If you want to use the FICO Community license and no longer want to see this message, use the following code before using the xpress module:
+ xpress.init('/home/al_dente/dev/xpress-demo/.venv/lib/python3.11/site-packages/xpress/license/community-xpauth.xpr')
+
+
+:: This is the solution
+FICO Xpress v9.2.5, Community, solve started 22:49:07, Nov 16, 2023
+Heap usage: 390KB (peak 390KB, 86KB system)
+Maximizing MILP noname using up to 4 threads and up to 7897MB memory, with these control settings:
+OUTPUTLOG = 1
+Original problem has:
+ 1 rows 16 cols 16 elements 16 entities
+Presolved problem has:
+ 1 rows 16 cols 16 elements 16 entities
+Presolve finished in 0 seconds
+Heap usage: 421KB (peak 433KB, 86KB system)
+
+Coefficient range original solved
+ Coefficients [min,max] : [ 2.67e+00, 4.89e+01] / [ 8.34e-02, 1.53e+00]
+ RHS and bounds [min,max] : [ 1.00e+00, 1.00e+02] / [ 1.00e+00, 3.12e+00]
+ Objective [min,max] : [ 1.06e+01, 2.67e+01] / [ 1.06e+01, 2.67e+01]
+Autoscaling applied standard scaling
+
+Will try to keep branch and bound tree memory usage below 7.0GB
+ *** Solution found: .000000 Time: 0.00 Heuristic: T ***
+ *** Solution found: 104.739705 Time: 0.00 Heuristic: e ***
+Starting concurrent solve with dual (1 thread)
+
+ Concurrent-Solve, 0s
+ Dual
+ objective dual inf
+ D 176.96135 .0000000
+------- optimal --------
+Concurrent statistics:
+ Dual: 1 simplex iterations, 0.00s
+Optimal solution found
+
+ Its Obj Value S Ninf Nneg Sum Dual Inf Time
+ 1 176.961353 D 0 0 .000000 0
+Dual solved problem
+ 1 simplex iterations in 0.00 seconds at time 0
+
+Final objective : 1.769613529197185e+02
+ Max primal violation (abs/rel) : 0.0 / 0.0
+ Max dual violation (abs/rel) : 0.0 / 0.0
+ Max complementarity viol. (abs/rel) : 0.0 / 0.0
+
+Starting root cutting & heuristics
+Deterministic mode with up to 1 additional thread
+
+ Its Type BestSoln BestBound Sols Add Del Gap GInf Time
+a 166.838660 176.961353 3 5.72% 0 0
+
+Performing root presolve...
+
+Reduced problem has: 1 rows 11 columns 11 elements
+Presolve dropped : 0 rows 5 columns 5 elements
+Will try to keep branch and bound tree memory usage below 6.9GB
+
+ Its Obj Value S Ninf Nneg Sum Dual Inf Time
+ 2 176.810476 D 0 0 .000000 0
+Optimal solution found
+Dual solved problem
+ 2 simplex iterations in 0.00 seconds at time 0
+
+Final objective : 1.768104758324006e+02
+ Max primal violation (abs/rel) : 0.0 / 0.0
+ Max dual violation (abs/rel) : 0.0 / 0.0
+ Max complementarity viol. (abs/rel) : 0.0 / 0.0
+
+Starting root cutting & heuristics
+Deterministic mode with up to 1 additional thread
+
+ Its Type BestSoln BestBound Sols Add Del Gap GInf Time
+ 1 K 166.838660 169.804364 3 1 0 1.75% 1 0
+
+Cuts in the matrix : 1
+Cut elements in the matrix : 11
+
+Performing root presolve...
+
+Reduced problem has: 2 rows 5 columns 10 elements
+Presolve dropped : 0 rows 6 columns 12 elements
+Presolve tightened : 5 elements
+Will try to keep branch and bound tree memory usage below 6.9GB
+
+ Its Obj Value S Ninf Nneg Sum Dual Inf Time
+ 4 170.072327 D 0 0 .000000 0
+Optimal solution found
+Dual solved problem
+ 4 simplex iterations in 0.00 seconds at time 0
+
+Final objective : 1.700723269425466e+02
+ Max primal violation (abs/rel) : 0.0 / 0.0
+ Max dual violation (abs/rel) : 0.0 / 0.0
+ Max complementarity viol. (abs/rel) : 0.0 / 0.0
+
+Starting root cutting & heuristics
+Deterministic mode with up to 1 additional thread
+
+ Its Type BestSoln BestBound Sols Add Del Gap GInf Time
+ *** Search completed ***
+Uncrunching matrix
+Final MIP objective : 1.668386599999999e+02
+Final MIP bound : 1.668388369613529e+02
+ Solution time / primaldual integral : 0.01s/ 24.267905%
+ Number of solutions found / nodes : 3 / 1
+ Max primal violation (abs/rel) : 0.0 / 0.0
+ Max integer violation (abs ) : 0.0
+166.83865999999995
+```
+
+
+
@@ 1,57 @@
+#!/usr/bin/python3
+
+import xpress as xp
+import pandas as pd
+import geopandas as gpd
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import numpy as np
+import folium
+
+def main():
+ gdf = gpd.read_file('knapsack_example.geojson')
+
+ print("\n:: This is the data")
+ print(gdf.shape)
+ print(gdf.head())
+
+ # decision variables from each Landscape Management Units (LMU)
+ x_i = [xp.var(i, vartype=xp.binary) for i in gdf['LMU_ID'].astype(str)]
+
+ # objective to maximize
+ Beta_i = gdf['res_d_SPM']
+ Z = xp.Sum(Beta_i * x_i)
+
+ # constraint (i.e. LMUs selected are constrained by number of acres)
+ alpha_i = gdf['Acres']
+ c1 = xp.Sum(alpha_i * x_i) <= 100 #try changing me to 80
+
+ # linear program
+ print("\n\n:: This is some licensing stuff, can be ignored")
+ p = xp.problem()
+ p.addVariable(x_i)
+ p.addConstraint(c1)
+ p.setObjective(Z, sense=xp.maximize)
+
+ print("\n\n:: This is the solution")
+ p.solve()
+ print(p.getObjVal())
+
+ # solution
+ var_solutions = [int(x.name) for x in x_i if p.getSolution(x) == 1]
+
+ # plot all LMUs and the objective variable
+ fig, ax = plt.subplots(figsize=(16,10))
+ gdf.plot(color='None', edgecolor='k', ax=ax)
+ gdf.plot(column='res_d_SPM', ax=ax, legend=True, scheme='quantiles', legend_kwds={'loc': 'lower right'})
+ ax.set_title('res_d_SPM', fontsize=18)
+
+ # highlight solution
+ outline_solutions = gdf.loc[gdf['LMU_ID'].isin(var_solutions)]
+ outline_solutions.plot(color='None', edgecolor='r', linewidth=5, ax=ax, hatch='///')
+
+ fig.savefig('choropleth_map.png')
+
+if __name__ == '__main__':
+ main()
+
