From b0614aefc663f0d454551d363eff3bf1c703e5df Mon Sep 17 00:00:00 2001 From: Andy Eschbacher Date: Mon, 8 May 2017 23:47:17 -0400 Subject: [PATCH] [mvp] fractional optimization --- src/pg/sql/25_optimization.sql | 18 +- .../crankshaft/analysis_data_provider.py | 40 +++++ .../crankshaft/optimization/optim.py | 160 ++++++++++-------- 3 files changed, 137 insertions(+), 81 deletions(-) diff --git a/src/pg/sql/25_optimization.sql b/src/pg/sql/25_optimization.sql index 42c892f..33f5012 100644 --- a/src/pg/sql/25_optimization.sql +++ b/src/pg/sql/25_optimization.sql @@ -1,11 +1,10 @@ -CREATE OR REPLACE FUNCTION -CDB_OptimAssignments(drain text, - source text, +CREATE OR REPLACE FUNCTION +CDB_OptimAssignments(source text, + drain text, drain_capacity text, source_production text, marginal_cost text, - waste_per_person numeric DEFAULT 0.01, - recycle_rate numeric DEFAULT 0.0, + dist_matrix_query text, dist_rate numeric DEFAULT 0.15, dist_threshold numeric DEFAULT null) RETURNS table(drain_id bigint, source_id int, cost numeric) AS $$ @@ -15,13 +14,12 @@ from crankshaft.optimization import Optim def cast_val(val): return float(val) if val is not None else None -params = {'waste_per_person': cast_val(waste_per_person), - 'recycle_rate': cast_val(recycle_rate), - 'dist_rate': cast_val(dist_rate), +params = {'dist_rate': cast_val(dist_rate), 'dist_threshold': cast_val(dist_threshold)} -optim = Optim(drain, source, drain_capacity, source_production, marginal_cost, - **params) + +optim = Optim(source, drain, dist_matrix_query, drain_capacity, + source_production, marginal_cost, **params) x = optim.output() return x diff --git a/src/py/crankshaft/crankshaft/analysis_data_provider.py b/src/py/crankshaft/crankshaft/analysis_data_provider.py index d733757..687e3e4 100644 --- a/src/py/crankshaft/crankshaft/analysis_data_provider.py +++ b/src/py/crankshaft/crankshaft/analysis_data_provider.py @@ -4,6 +4,8 @@ import pysal_utils as pu import numpy as np class AnalysisDataProvider(object): + """Analysis providers for crankshaft functions. These rely on database + access through `plpy`""" def get_getis(self, w_type, params): """fetch data for getis ord's g""" try: @@ -91,6 +93,44 @@ class AnalysisDataProvider(object): resp = plpy.execute(query) return np.array(resp[0]['col'], dtype=dtype) + def get_distance_matrix(self, query, origin_ids, destination_ids): + """Transforms a SQL table origin-destination table into a distance + matrix. + + :param query: Query that exposes the data needed for building the + distance matrix. Query should have the following columns: + - origin_id (int) + - destination_id (int) + - length_km (numeric) + :type query: str + :param origin_ids: List of origin IDs + :type origin_ids: list of ints + :param destination_ids: List of origin IDs + :type destination_ids: list of ints + :returns: 2D array of distances from all origins to all destinations + :rtype: numpy.array + """ + try: + resp = plpy.execute(''' + SELECT "origin_id", "destination_id", "length_km" + FROM ({query}) as _wrap + '''.format(query=query)) + except plpy.SPIError as err: + plpy.error("Failed to build distance matrix: {}".format(err)) + + pairs = {(row['origin_id'], row['destination_id']): row['length_km'] + for row in resp} + distance_matrix = np.array([ + pairs[(origin, destination)] + for destination in destination_ids + for origin in origin_ids + ]) + + return np.array(distance_matrix, + dtype=float).reshape((len(destination_ids), + len(origin_ids))) + + def get_pairwise_distances(self, drain_query, source_query, id_col='cartodb_id'): """Retuns the pairwise distances between row i and j for all i in diff --git a/src/py/crankshaft/crankshaft/optimization/optim.py b/src/py/crankshaft/crankshaft/optimization/optim.py index ba91b69..e2d0639 100644 --- a/src/py/crankshaft/crankshaft/optimization/optim.py +++ b/src/py/crankshaft/crankshaft/optimization/optim.py @@ -9,47 +9,44 @@ class Optim(object): """Linear optimization class for logistics cost minimization Optimization for logistics based on models: - - amount_per_unit * (1 - recycle_rate) * population - source_amount * (marginal_cost + transport_cost * distance) - That is, `cost ~ population * distance` """ - def __init__(self, drain_query, source_table, capacity_column, # pylint: disable=too-many-arguments - production_column, marginal_column, **kwargs): + def __init__(self, source_query, drain_query, dist_matrix_query, + capacity_column, production_column, marginal_column, + **kwargs): # set data provider - defaults to SQL database access self.data_provider = kwargs.get('data_provider', AnalysisDataProvider()) # model parameters self.model_params = { - 'amount_per_unit': kwargs.get('amount_per_unit', 0.01), 'dist_cost': kwargs.get('dist_cost', 0.15), - 'recycle_rate': kwargs.get('recycle_rate', 0.0), - 'dist_threshold': kwargs.get('dist_threshold', None)} + 'dist_threshold': kwargs.get('dist_threshold', None), + 'solver': kwargs.get('solver', 'glpk')} self._check_model_params() - # model data - self.model_data = { - 'drain_capacity': self.data_provider.get_column(drain_query, - capacity_column), - 'source_amount': (self.model_params['amount_per_unit'] * - (1. - self.model_params['recycle_rate']) * - self.data_provider.get_column(source_table, - production_column)), - 'marginal_cost': self.data_provider.get_column(drain_query, - marginal_column), - 'distance': self.data_provider.get_pairwise_distances(source_table, - drain_query) - } - self.model_data['cost'] = self.calc_cost() # database ids self.ids = { 'drain': self.data_provider.get_column(drain_query, 'cartodb_id', dtype=int), - 'source': self.data_provider.get_column(source_table, + 'source': self.data_provider.get_column(source_query, 'cartodb_id', dtype=int)} + # model data + self.model_data = { + 'drain_capacity': self.data_provider.get_column(drain_query, + capacity_column), + 'source_amount': self.data_provider.get_column(source_query, + production_column), + 'marginal_cost': self.data_provider.get_column(drain_query, + marginal_column), + 'distance': + self.data_provider.get_distance_matrix(dist_matrix_query, + self.ids['source'], + self.ids['drain'])} + self.model_data['cost'] = self.calc_cost() self.n_sources = len(self.ids['source']) self.n_drains = len(self.ids['drain']) @@ -73,22 +70,17 @@ class Optim(object): def _check_model_params(self): """Ensure model parameters are well formed""" - if (self.model_params['recycle_rate'] is None or - self.model_params['recycle_rate'] < 0 or - self.model_params['recycle_rate'] > 1): - raise ValueError("`recycle_rate` must be between 0 and 1.") - - if (self.model_params['amount_per_unit'] is None or - self.model_params['amount_per_unit'] < 0): - raise ValueError("`amount_per_unit` must be greater than zero.") - if (self.model_params['dist_threshold'] <= 0 and self.model_params['dist_threshold'] is not None): raise ValueError("`dist_threshold` must be greater than zero") if (self.model_params['dist_cost'] is None or self.model_params['dist_cost'] < 0): - raise ValueError("`dist_cost must be greater than zero") + raise ValueError("`dist_cost` must be greater than zero") + + if self.model_params['solver'] not in (None, 'glpk'): + raise ValueError("`solver` must be one of 'glpk' (default) " + "or None.") return None @@ -99,8 +91,7 @@ class Optim(object): transport from source to drain :rtype: List of tuples """ - - # n_drains x n_sources matrix (row, column) + # retrieve fractional assignments assignments = self.optim() # crosswalks for matrix index -> cartodb_id @@ -118,11 +109,11 @@ class Optim(object): nonzeros = np.nonzero(assignments) source_index, drain_index = nonzeros[0], nonzeros[1] # - assigned_costs = [(drain_id_crosswalk[drain_index[i]], - source_id_crosswalk[source_index[i]], - self.model_data['cost'][drain_index[i], - source_index[i]]) - for i in range(len(source_index))] + assigned_costs = [(drain_id_crosswalk[drain_index[idx]], + source_id_crosswalk[source_val], + self.model_data['cost'][drain_index[idx], + source_val]) + for idx, source_val in enumerate(source_index)] return assigned_costs def cost_func(self, distance, waste, marginal): @@ -148,12 +139,13 @@ class Optim(object): Populate an d x s matrix according to the cost equation :returns: d x s matrix of costs from area i to plant j - :rtype: NumPy matrix + :rtype: numpy.array """ - costs = np.array([self.cost_func(distance, - self.model_data['source_amount'][pair[1]], - self.model_data['marginal_cost'][pair[0]]) - for pair, distance in np.ndenumerate(self.model_data['distance'])]) + costs = np.array( + [self.cost_func(distance, + self.model_data['source_amount'][pair[1]], + self.model_data['marginal_cost'][pair[0]]) + for pair, distance in np.ndenumerate(self.model_data['distance'])]) return costs.reshape(self.model_data['distance'].shape) def optim(self): @@ -163,51 +155,77 @@ class Optim(object): minimize c'*x by assigning x values subject to G*x <= h A*x = b - x[k] is binary - :returns: Assignments array (of 1s and 0s) of shape c.T - :rtype: NumPy array + 0 <= x[k] <= 1 + :returns: Fractional assignments array (of 1s and 0s) of shape c.T. + Value at position (i, j) corresponds to the fraction of source + `i`'s supply to drain `j`. + + :rtype: numpy.array """ + n_pairings = self.n_sources * self.n_drains + # --- # costs # elements chosen to minimize sum - cost = cvxopt.matrix(self.model_data['cost'].ravel('F')) + cost = np.nan_to_num(self.model_data['cost']) + cost = cvxopt.matrix(cost.ravel('F')) # --- # equality constraint variables # each area is serviced once - A = cvxopt.spmatrix(1., # pylint: disable=invalid-name + A = cvxopt.spmatrix(1., [i // self.n_drains - for i in range(self.n_drains * self.n_sources)], - range(self.n_drains * self.n_sources)) - b = cvxopt.matrix(np.ones((self.n_sources, 1)), tc='d') # pylint: disable=invalid-name + for i in range(n_pairings)], + range(n_pairings), tc='d') + b = cvxopt.matrix([1.] * self.n_sources, tc='d') + + # make nan's in cost impossible + if np.isnan(self.model_data['distance']).any(): + i_vals, j_vals = np.where(np.isnan(self.model_data['distance'])) + for idx, i_val in enumerate(i_vals): + i = int(i_val) + j = int(i_val * self.n_drains + j_vals[idx]) + A[i, j] = 0 # knock out values above distance threshold if self.model_params['dist_threshold']: - j_locs, i_locs = np.where(self.model_data['distance'] > 100) - for idx, ival in enumerate(i_locs): - A[int(ival), int(ival * 10 + j_locs[idx])] = 0 + j_vals, i_vals = np.where(self.model_data['distance'] > + self.model_params['dist_threshold']) + for idx, ival in enumerate(i_vals): + A[int(ival), int(ival * self.n_drains + j_vals[idx])] = 0 # --- # inequality constraint variables # each plant never goes over capacity - drain_capacity = cvxopt.matrix(self.model_data['drain_capacity'], - tc='d') - source_amounts = cvxopt.spmatrix( - np.repeat(self.model_data['source_amount'], self.n_drains), - [i % self.n_drains for i in range(self.n_drains * self.n_sources)], - range(self.n_drains * self.n_sources)) + drain_capacity = cvxopt.matrix([ + cvxopt.matrix(self.model_data['drain_capacity'], tc='d'), + cvxopt.matrix([1.] * n_pairings, tc='d'), + cvxopt.matrix([0.] * n_pairings, tc='d') + ]) - binary_entries = set(range(self.n_drains * self.n_sources)) + # inequality maxima + ineq_maxs = cvxopt.sparse([ + cvxopt.spmatrix( + np.repeat(self.model_data['source_amount'], self.n_drains), + [i % self.n_drains for i in range(n_pairings)], + range(n_pairings), tc='d'), + cvxopt.spmatrix(1., + range(n_pairings), + range(n_pairings)), + cvxopt.spmatrix(-1., + range(n_pairings), + range(n_pairings)) + ], tc='d') # solve - (sol, assignments) = ilp(c=cost, G=source_amounts, h=drain_capacity, - A=A, b=b, B=binary_entries) - if sol != 'optimal': + sol = cvxopt.solvers.lp(c=cost, G=ineq_maxs, h=drain_capacity, + A=A, b=b, solver=self.model_params['solver']) + if sol['status'] != 'optimal': raise Exception("No solution possible: {}".format(sol)) - assign_shape = (self.model_data['cost'].shape[1], - self.model_data['cost'].shape[0]) - - # Note: assignments needs to be shaped like self.model_data['cost'].T - return np.array(assignments, - dtype=int).flatten().reshape(assign_shape) + # NOTE: assignments needs to be shaped like self.model_data['cost'].T + return np.array(sol['x'], + dtype=float)\ + .flatten()\ + .reshape((self.model_data['cost'].shape[1], + self.model_data['cost'].shape[0]))