Uso de la programación lineal para la resolución de un modelo de flujos de cuota de mercado

Uso de la programación lineal para la resolución de un modelo de flujos de cuota de mercado

Francisco Espiga Fernández
Enero, 2022

Más sobre ciencia de datos: cienciadedatos.net

En este artículo exploraremos el uso de modelos de programación lineal sobre los datos de partida de la cuota de mercado de un conjunto de competidores, para determinar los flujos individuales de captación de cuota entre ellos.

Caso de negocio


En muchos contextos e industrias nos encontramos con un conjunto de productos o empresas que compiten entre ellas por la cuota de mercado. Ejemplos de esto podemos encontrarlos en:

  • La cuota por provincia en Estaciones de Servicio en la industria de Oil&Gas, donde CORES en España es una buena aproximación de ese share.

  • Las audiencias de las distintas cadenas de televisión.

  • El share of wallet de distintas marcas del consumidor y que recopilan empresas como Kantar o Statista.

Algo común en todos estos casos, es que, si bien conocemos el valor preciso o una estimación de estos market share, desconocemos el trasvase de cuota entre los distintos competidores, para cada instante de tiempo, así como las dinámicas que favorecen o dificultan dicho trasvase. En este artículo, recurriremos a los grafos como herramienta de visualización y a un modelo de programación lineal, o LP (Linear Programming).

Antes de seguir avanzando, debemos indicar al lector que para poder ejecutar correctamente el código del artículo, es necesario instalar la librería pyomo así como el solver cbc del proyecto COIN-OR

Para la librería, se puede instalar tanto via pip como desde conda independientemente de nuestro sistema operativo. Para el solver, dependiendo del sistema operativo, recomendamos:

  • WINDOWS: descargar de la web del proyecto COIN-OR el solver, e indicar la ruta del ejecutable pyo.SolverFactory('cbc', executable = 'path/to/solver/cbc.exe')

  • LINUX: instalación via conda mediante el comando conda install -c conda-forge coincbc

Definiendo el problema


En el caso que nos ocupa, vamos a utilizar un modelo de programación lineal, que modelaremos usando la librería pyomo para determinar los flujos de trasvase de cuotas de mercado entre competidores y recurriremos a la librería networkx para la visualización de los datos en forma de grafo.

Tomaremos para ello los datos de cuota de mercado a lo largo de un año, para los competidores. Dado que nuestro caso de uso podría ser de aplicación en numerosos ámbitos, generaremos nuestro escenario de datos de manera aleatoria y para un número arbitrariamente escogido de competidores de 10.

Generación y visualización de los datos


Generaremos los datos de cuota de mercado para todo el año 2020 y para los 10 competidores que hemos escogido de manera arbitraria. Como lo que queremos calcular son los flujos de entrada de cuota de mercado, calcularemos también con pandas el diferencial mes a mes de cuota de mercado para cada competidor.

In [1]:
# Librerías
# ======================================================================================
import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt
from tqdm import tqdm 
import networkx as nx
import pyomo.environ as pyo
from pyomo.core.base import Constraint
from pyomo.opt import SolverStatus, TerminationCondition
In [2]:
#!conda install -c conda-forge coincbc
In [3]:
# Generación de un escenario aleatorio
# ======================================================================================
n_competitors = 10
start_date = '2020-01-01'
end_date = '2021-01-01'
print(f"Generating random scenario from {start_date} to {end_date} with {n_competitors} competitors")

dates = [i.strftime('%Y-%m-%d') for i in pd.date_range(start_date, end_date, freq = 'M')]
CPs = [f"CP_{i}" for i in range(8)]

data = pd.DataFrame(np.random.uniform(0,100, [len(dates), len(CPs)]), index = dates, columns = CPs)
data_mktshare = data.div(data.sum(axis = 1).values, axis = 0)
data_delta_mktshare = data_mktshare.diff().fillna(value=0)
Generating random scenario from 2020-01-01 to 2021-01-01 with 10 competitors

Vamos a visualizar ahora la cuota de mercado de los competidores en un grafo, para un mes al azar.

In [4]:
# Gráfico del grafo
# ======================================================================================
timestamp = np.random.choice(data_mktshare.index, 1)[0]
attr = 'mkt_share'        
edgelist = pd.DataFrame({'source':CPs}).merge(pd.DataFrame({'target':CPs}), how = 'cross')
G = nx.from_pandas_edgelist(edgelist, source='source', target='target')
        
nx.set_node_attributes(G, data_mktshare[data_mktshare.index == timestamp].reset_index(drop=True).T.to_dict()[0], attr)


COLOR_SCHEME = "RdBu" 
colors = [G.nodes[node][attr] for node in list(G.nodes())]
nx.draw_networkx(G, node_color=colors, cmap=COLOR_SCHEME)#, vmin = -1, vmax = 1)

ax=plt.gca()
PCM=ax.get_children()[0]  
plt.colorbar(PCM, ax=ax)
plt.show()

Verificamos como algunos de los competidores tienen un porcentaje muy elevado de cuota de mercado mientras que otros están un orden de magnitud por debajo, así como que todos nuestros datos son positivos y suman uno.

Modelización

En primer lugar, tenemos que hacer una serie de asunciones y simplificaciones para nuestro problema.

Asunciones

  • El flujo entrante de cuota de mercado de un competidor, proviene del resto de competidores. En nuestro contexto parece trivial, pero por ejemplo en el financiero, podríamos intentar ver el flujo monetario entre distintos competidores, con una fuente de financiación externa en cash.

  • Denominamos flujos de entrada o inflows a los flujos positivos, y de salida o outflows a los negativos.

  • No puede haber bucles, es decir, flujos de entrada o salida de un competidor hacia sí mismo.

Sets

Comenzamos por definir los distintos conjuntos de nuestro problema. En el contexto de la modelización y programación lineal, son los índices en los que indexaremos las distintas variables, parámetros y restricciones. En nuestro caso, tenemos los competidores $C$ y los meses $M$.

  • $C$ : los competidores del mercado.

  • $T$: los meses de histórico.

Parámetros

Se denominan parámetros a aquellos valores que nos son conocidos de antemano y sobre cuyos valores no decidimos durante la optimización. En nuestro caso, tanto el marketshare de un competidor en un mes, como la evolución mes a mes del mismo.

  • $\Delta Mk_{c,m}$: la evolución del marketshare entre el mes $m$ y el mes $m+1$ de un competidor $c$.

  • $Mk_{c,m}$: el marketshare de un competidor en un mes dado.

Variables de decisión

En este caso, nuestras variables de decisión van a ser los flujos de entrada o salida de cuota de mercado entre competidores. Podríamos adoptar la convención de un flujo único que pudiera ser tanto positivo ($Flow_{a\rightarrow b}$) si $b$ captura cuota de mercado a $a$, como negativo en caso contrario. Sin embargo, hemos optado por desdoblarlos para simplificar la redacción de algunas de las restricciones del modelo, y dotarlo de mayor flexibilidad en caso de querer expandirlo en un futuro.

  • $\forall m \in T, c_1 \in C, c_2 \in C: OutFlow_{m,c_1,c_2} \in R^+$: flujo negativo de cuota de mercado entre competidores.

  • $\forall m \in T, c_1 \in C, c_2 \in C: InFlow_{m,c_1,c_2} \in R^+$: flujo positivo de cuota de mercado entre competidores.

Consideraciones adicionales

Formato de los datos de entrada

Para inicializar los valores de los parámetros en pyomo, necesitamos tener los datos en formato largo, con las columnas índice (competidores y meses) y la columna valor. Esto lo usaremos para declarar un diccionario con claves las tuplas de los índices y valor el valor del parámetro que queremos inicializar, o bien filtrar la fila deseada del dataframe.

In [5]:
data_long_delta_mktshare = (
    data_delta_mktshare
    .reset_index()
    .melt(id_vars='index')
    .rename(columns={'variable':'CP', 'value':'delta_mkt_share', 'index':'timestamp'})
)

data_long_mktshare = (
    data_mktshare
    .reset_index()
    .melt(id_vars='index')
    .rename(columns={'variable':'CP', 'value':'mkt_share', 'index':'timestamp'})
)

data_long_mktshare.head()
Out[5]:
timestamp CP mkt_share
0 2020-01-31 CP_0 0.246785
1 2020-02-29 CP_0 0.191869
2 2020-03-31 CP_0 0.131337
3 2020-04-30 CP_0 0.029125
4 2020-05-31 CP_0 0.085638

Variables de holgura

Añadimos dos variables de holguras, vPositiveSlack y vNegativeSlack para facilitar el determinar los flujos. En nuestro caso de programación lineal, puede que algunas restricciones, como la de simetría (el flujo de salida de $B$ a $A$ es igual al de entrada desde $B$ a $A$, cambiado de signo) puede que no sean exactamente iguales, sino con tolerancias cercanas al cero. En caso de no ser idénticos, esto haría que nuestro problema no encontrase una solución que respetase todas las restricciones, es decir, que fuese factible.

Al recurrir a las holguras, damos cierto "margen" para poder cumplir las restricciones, y dotamos al problema de flexibilidad en caso de expandirlo en el futuro con relaciones y dependencias más complejas, pero garantizando que exista una solución factible. Sin embargo, es necesario que penalicemos el uso de estas variables de holgura, ya que si no, nuestro problema no tendría ningún incentivo en respetar las restricciones. Esto lo lograremos en la función objetivo.

Restricciones implícitas en el procesamiento de datos

En nuestros datos generados, así como en un contexto que utilice datos reales, existen dos restricciones de manera implícita que no necesitamos indicar en el modelo. Por un lado, que la suma de las cuotas de mercado entre todos los competidores sea del 100% (sume 1) y, por otro lado, que la suma de las evoluciones de las cuotas de mercado ($\Delta Mk$) sea nulo.

  • CTa: Cálculo del delta marketshare: $$\forall c\in C, m\in T-\{M_0\}:\Delta Mk_{z,t} = Mk_{c,m}-Mk_{c, m-1}$$

  • CTb: Ecuación de equilibrio de la cuotrra de mercado: $$\forall m\in T: \sum_c Mk_{c,m} == 1$$

In [6]:
# Optimización del modelo
# ======================================================================================

model = pyo.ConcreteModel()

# Sets ----
model.sCP = pyo.Set(initialize = CPs)
model.sTimestamps = pyo.Set(initialize = dates)

# Parameters ----
def get_mktshare_delta(model, p, m):
    return data_long_delta_mktshare[(data_long_delta_mktshare.CP == p) & (data_long_delta_mktshare.timestamp == m)]['delta_mkt_share'].values[0]

model.pDeltaMktShare = pyo.Param(model.sCP, model.sTimestamps, initialize = get_mktshare_delta)

def get_mktshare(model, p, m):
    return data_long_mktshare[(data_long_mktshare.CP == p) & (data_long_mktshare.timestamp == m)]['mkt_share'].values[0]
    
model.pMktShare = pyo.Param(model.sCP, model.sTimestamps, initialize = get_mktshare)

# Variables ----
"""
Outflow A,B es lo que el nodo B transfiere al nodo A.
Inflow A,B es lo que el nodo A toma del nodo B. 
Por tanto: Outflow A,B = -Inflow B,A

Para un nodo: la variación de marketshare es la suma de los outflows (A, x): todo lo que se toma de los nodos; y los inflows (A, y) que es recibido por parte de los demás nodos.
"""
model.vOutFlow = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.NonPositiveReals) # Negativo para flujos de salida.
model.vInFlow = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.NonNegativeReals) # Positivo para flujos entrantes.
model.vPositiveSlack = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.NonNegativeReals, bounds=(0,0.05))
model.vNegativeSlack = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.NonPositiveReals, bounds=(-0.05,0))

Restricciones

CT1: Market delta flows

La variación del marketshare en un mes dado es la suma de los flujos de entrada y de salida respecto del resto de competidores.

$$\forall c\in C, t\in T:\Delta M_{c,t} = \sum_{y\in C-\{c\}} OutFlow_{t,y,c} + InFlow_{t,y,c} $$

In [7]:
## CT1: el delta de marketshare debe ser la suma de inflows (>=0) y outflows (<=0)
## Añadimos el slack u holgura para facilitar encontrar una solución factible. Nuestro objetivo será minimizarlo.
def mktshare_equilibrium(model, p, t):
    return model.pDeltaMktShare[p, t]==sum(model.vOutFlow[p, p_x, t] + model.vInFlow[p, p_x, t] + model.vPositiveSlack[p,p_x, t] + model.vNegativeSlack[p,p_x, t] for p_x in model.sCP)

model.ctMktShareEquilibrium = pyo.Constraint(model.sCP, model.sTimestamps, rule = mktshare_equilibrium)

CT2: Out/In flow balance:

El flujo de entrada de A desde B, es el negativo del flujo de salida de B hacia A. $$\forall t \in T, z_1 \in Z, z_2 \in Z: OutFlow_{t,z_1,z_2}==InFlow_{t,z_2,z_1}$$

In [8]:
## CT2: Simetría Outflow (A,B) e Inflow (B,A), cambiando el signo. 
def connect_flows(model, p_a, p_b, t):
    return model.vOutFlow[p_a, p_b, t]==-model.vInFlow[p_b, p_a, t]
model.ctFlowConnection = pyo.Constraint(model.sCP, model.sCP, model.sTimestamps, rule = connect_flows)

CT3: Evitar bucles

  • CT3a: salida: $$\forall t \in T, z_1 \in Z, z_1 \in Z: OutFlow_{t,z_1,z_1}=0$$
  • CT3b: entrada: $$\forall t \in T, z_1 \in Z, z_1 \in Z: InFlow_{t,z_1,z_1}=0$$
In [9]:
## CT3: Evitar bucles de entrada y de salida
def outflow_loops(model, p, t):
    return model.vOutFlow[p, p, t]==0

model.ctBlockOutflowLoops = pyo.Constraint(model.sCP, model.sTimestamps, rule = outflow_loops)

def inflow_loops(model, p, t):
    return model.vInFlow[p, p, t]==0
    
model.ctBlockInflowLoops = pyo.Constraint(model.sCP, model.sTimestamps, rule = inflow_loops)

CT4: limitar los flujos de salida

Los flujos de salida de un competidor hacia el resto, siempre deben ser menores que el marketshare del mes previo. $$\forall t \in T, z_1 \in Z, z_2 \in Z: OutFlow_{t,z_1,z_2}<=MktShare_{t-1, z_1}$$

In [10]:
## CT4: La suma de flujos de salida debe ser siempre inferior al marketshare precedente.
def sum_of_outflows(model, p, t):
    if t == model.sTimestamps.first():
        return Constraint.Skip
    else:
        return -(sum(model.vOutFlow[p, p_x, t] for p_x in model.sCP))<=model.pMktShare[p, model.sTimestamps.prev(t)]
        
model.ctMaxOutflows = pyo.Constraint(model.sCP, model.sTimestamps, rule = sum_of_outflows)

CT5: unicidad de los flujos de entrada y de salida

En nuestro contexto de negocio, no tendría sentido que un competidor A capture cuota de mercado de un competidor B, y que el competidor B capture cuota de mercado del competidor A. Podríamos argumentar en contra, que si por ejemplo el trasvase de cuota de mercado obedece a criterios económicos por un lado; y por otro a la presencia en RRSS, publicidad en medios, etc, podría darse el caso de que un competidor A captura "económicamente" cuota de B, y éste "por publicidad", cuota de A. No obstante, en nuestro caso mantenemos un enfoque simple de flujos de entrada o de salida en términos netos, por lo que no se puede dar esta situación.

Recordemos que los flujos de entrada o inflows tienen el signo positivo, y los outflows o de salida, negativos. En este caso, podemos restringir esa unicidad de los flujos del siguiente modo:

  • Para los inflows: $Inflow_{a,b}+Inflow_{b,a}=max(Inflow_{a,b},Inflow_{b,a})$
  • Para los outflows: $Outflow_{a,b}+Outflow_{b,a}=min(Outflow_{a,b},Outflow_{b,a})$

No obstante, tanto el máximo como el mínimo no son funciones lineales, pero las podemos modelar mediante el uso de variables binarias y restricciones de tipo "big M". Explicamos el caso de los outflows en detalle, aunque para los inflows será análogo.

  • El mínimo siempre será $\leq$ que ambos flujos de salida, al ser el más negativo de los dos.
$$Minimum\;outflow_{t, a, b}= Min(Outflow_{t,a,b}, Outflow_{t,b,a})$$$$Minimum\;outflow_{t, a, b}<= Outflow_{t,a,b}$$$$Minimum\;outflow_{t, a, b}<= Outflow_{t,b,a}$$

  • Por otro lado reforzamos que la suma de ambos tenga que ser igual al mínimo, de manera simétrica.
$$Minimum\;outflow_{t, a, b}== Outflow_{t,a,b} + Outflow_{t,b,a}$$

$$Minimum\;outflow_{t, a, b}==Minimum\;outflow_{t, b, a}$$

Para simplificar la notación, nos referimos a los flujos de salida como x1 y x2; al valor mínimo como X, M como un número real positivo y un orden de magnitud mayor que nuestros valores (e.g. 20) e y a una variable binaria o indicativa. Tenemos:

$$X<=x1$$$$X<=x2$$$$X>=x1-M(1-y)$$$$X>=x2-My$$

Las dos primeras restricciones son las anteriormente mencionadas. Las dos siguientes, con el uso de la constante M y la variable binaria indicativa y, fuerzan que solo uno de los valores pueda tomar un valor distinto de cero.

El planteamiento es análogo para los flujos de salida, reemplazando el mínimo por el máximo y las convenciones de signos de las inecuaciones.

In [11]:
M = 1e2

## CT5: bien hay un flujo de salida de A a B o bien de B a A.
"""
No podemos tener flujos de salida o de entrada de A a B y de B a A simultáneamente. O bien uno transfiere (recibe) o bien el otro. 
Esto es lo mismo matemáticamente que decir:
- La suma de outflows es igual al outflow minimo.
- La sunma de inflows es igual al inflow maximo.

    x2-x1 = My
    x1-x2 = M(1-y)
    X<=x1
    X<=x2
    X>=x1-M(1-y)
    X>=x2-My
"""
### El mínimo outflow es la suma de los outflows.
### ref: https://or.stackexchange.com/questions/1160/how-to-linearize-min-function-as-a-constraint
model.vMinOutflow = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.NonPositiveReals)
model.vbMinOutflowFlag = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.Binary)

model.ctOneOutflow = pyo.ConstraintList()
for p_a in model.sCP:
    for p_b in model.sCP:
        for t in model.sTimestamps:
            #model.ctOneOutflow.add(expr = model.vOutFlow[p_b, p_a, t]-model.vOutFlow[p_a, p_b, t] == M*model.vbMinOutflowFlag[p_a, p_b, t])
            #model.ctOneOutflow.add(expr = model.vOutFlow[p_a, p_b, t]-model.vOutFlow[p_b, p_a, t] == M*(1-model.vbMinOutflowFlag[p_a, p_b, t]))

            model.ctOneOutflow.add(expr = model.vMinOutflow[p_a, p_b, t] <= model.vOutFlow[p_a, p_b, t])
            model.ctOneOutflow.add(expr = model.vMinOutflow[p_a, p_b, t] <= model.vOutFlow[p_b, p_a, t])

            model.ctOneOutflow.add(expr = model.vMinOutflow[p_a, p_b, t] >= model.vOutFlow[p_a, p_b, t]-M*(1-model.vbMinOutflowFlag[p_a, p_b, t]))
            model.ctOneOutflow.add(expr = model.vMinOutflow[p_a, p_b, t] >= model.vOutFlow[p_b, p_a, t]-M*(model.vbMinOutflowFlag[p_a, p_b, t]))

            model.ctOneOutflow.add(expr = model.vOutFlow[p_b, p_a, t]+model.vOutFlow[p_a, p_b, t] == model.vMinOutflow[p_a, p_b, t])
            model.ctOneOutflow.add(expr = model.vMinOutflow[p_b, p_a, t] == model.vMinOutflow[p_a, p_b, t])


### El máximo outflow es la suma de los outflows.
### ref: https://math.stackexchange.com/questions/2446606/linear-programming-set-a-variable-the-max-between-two-another-variables
model.vMaxInflow = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.NonNegativeReals)
model.vbMaxInflowFlag = pyo.Var(model.sCP, model.sCP, model.sTimestamps, within = pyo.Binary)

model.ctOneInflow = pyo.ConstraintList()
for p_a in model.sCP:
    for p_b in model.sCP:
        for t in model.sTimestamps:
            model.ctOneInflow.add(expr = model.vMaxInflow[p_a, p_b, t] >= model.vInFlow[p_a, p_b, t])
            model.ctOneInflow.add(expr = model.vMaxInflow[p_a, p_b, t] >= model.vInFlow[p_b, p_a, t])

            model.ctOneInflow.add(expr = model.vMaxInflow[p_a, p_b, t] <= model.vInFlow[p_a, p_b, t] + M*model.vbMaxInflowFlag[p_a, p_b, t])
            model.ctOneInflow.add(expr = model.vMaxInflow[p_a, p_b, t] <= model.vInFlow[p_b, p_a, t] + M*(1-model.vbMaxInflowFlag[p_a, p_b, t]))

            model.ctOneInflow.add(expr = model.vInFlow[p_b, p_a, t]+model.vInFlow[p_a, p_b, t] == model.vMaxInflow[p_a, p_b, t])
            model.ctOneInflow.add(expr = model.vMaxInflow[p_b, p_a, t] == model.vMaxInflow[p_a, p_b, t])

CT6: restricción de inicialización de los flujos en T_0

Como última restricción, obligamos a que los flujos de entrada y salida para todos los competidores en el instante 0, deban ser nulos. Esto garantizará también el cumplimiento de las restricciones de máximo / mínimo anteriormente formuladas

$$\forall c_1 \in C, c_2 \in C: InFlow_{\;0,c_1,c_2}=0$$$$\forall c_1 \in C, c_2 \in C: OutFlow_{\;0,c_1,c_2}=0$$
In [12]:
## CT6: Los flujos en T=0 son nulos.
def zero_inflows(model, p_a, p_b):
    return model.vInFlow[p_a, p_b, model.sTimestamps.first()] ==0
model.ctZeroInitialInFlows = pyo.Constraint(model.sCP, model.sCP, rule = zero_inflows)

def zero_outflows(model, p_a, p_b):
    return model.vOutFlow[p_a, p_b, model.sTimestamps.first()] ==0
model.ctZeroInitialOutFlows = pyo.Constraint(model.sCP, model.sCP, rule = zero_outflows)

Función objetivo

Por último, nuestra función objetivo no es trivial. Aunque el caso que nos ocupa, al querer repercutir las variaciones de cuota de mercado a flujos de entrada y salida a competidores, podría formularse como un problema de satisfacción de restricciones, nosotros hemos añadido las variables de holgura o slacks para facilitar el encontrar una solución factible y tener flexibilidad de cara a modificaciones ulteriores. Por tanto, nuestra función objetivo será la suma de las holguras positivas, menos la suma de las holguras negativas, que buscaremos minimizar.

$$min z = \sum_{t, z_1, z_2} slack^+_{t, z_1,z_2} - \sum_{t, z_1, z_2} slack^-_{t, z_1,z_2}$$
In [13]:
# Función objetivo
minimize = 1
model.Z = pyo.Objective(sense = minimize, expr = sum(
    model.vPositiveSlack[p_a,p_b,t] for p_a in model.sCP for p_b in model.sCP for t in model.sTimestamps)-sum(
    model.vNegativeSlack[p_a,p_b,t] for p_a in model.sCP for p_b in model.sCP for t in model.sTimestamps))

Solucionando el modelo e inspeccionando los resultados


Una vez concluida la fase de modelado, podemos proceder a resolverlo y analizar y visualizar los resultados. En este caso, usaremos el solver cbc del proyecto COIN-OR, un solver lineal y open source. Como alternativas gratuitas a éste, SCIP recientemente ha sido disponibilizado tanto para entornos académicos como para profesionales sin necesidad de licenciamiento.

In [14]:
solver = pyo.SolverFactory('cbc')

solution_status =  solver.solve(model)
solution_status
Out[14]:
{'Problem': [{'Name': 'unknown', 'Lower bound': 0.0, 'Upper bound': 0.0, 'Number of objectives': 1, 'Number of constraints': 4489, 'Number of variables': 3314, 'Number of nonzeros': 850, 'Sense': 'minimize'}], 'Solver': [{'Status': 'ok', 'User time': -1.0, 'System time': 1.13, 'Wallclock time': 1.13, 'Termination condition': 'optimal', 'Termination message': 'Model was solved to optimality (subject to tolerances), and an optimal solution is available.', 'Statistics': {'Branch and bound': {'Number of bounded subproblems': 0, 'Number of created subproblems': 0}, 'Black box': {'Number of iterations': 0}}, 'Error rc': 0, 'Time': 1.5465221405029297}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

El modelo es capaz de encontrar rápidamente una solución factible y óptima. Uno de los mayores desafíos en el ámbito de la programación lineal, y que hace que en entornos industriales debamos decantarnos a favor de solvers comerciales, es su capacidad para encontrar rápidamente (en términos relativos) soluciones que sean factibles, sobre las que poder continuar la exploración de cara a encontrar una óptima. En las ocasiones en que utilicemos un solver open source, y experimentemos largos tiempos de resolución que no logran encontrar una solución, o cuyo resultado sea el de infactibilidad, puede resultar conveniente el uso de dichos solvers comerciales, ya que ellos sí que pueden ser capaces de encontrarla, o de arrojar el resultado de la infactibilidad con mayor rapidez.

In [15]:
def prepare_inflows(model, exclude_zeros:bool = True):
    inflows = pd.DataFrame.from_dict(model.vInFlow.extract_values(), orient = 'index').reset_index()
    if exclude_zeros:
        inflows = inflows[inflows[0]>0]

    inflows[['cp_a', 'cp_b', 'timestamp']] = pd.DataFrame(inflows['index'].tolist(), index=inflows.index)
    inflows = inflows.drop('index', axis = 1).rename(columns = {0:'inflow', 'cp_a':'target', 'cp_b':'source'})
    return inflows

def prepare_outflows(model, exclude_zeros:bool = True):
    outflows = pd.DataFrame.from_dict(model.vOutFlow.extract_values(), orient = 'index').reset_index()
    if exclude_zeros:
        outflows = outflows[outflows[0]<0]

    outflows[['cp_a', 'cp_b', 'timestamp']] = pd.DataFrame(outflows['index'].tolist(), index=outflows.index)
    outflows = outflows.drop('index', axis = 1).rename(columns = {0:'outflow', 'cp_a':'source', 'cp_b':'target'})
    return outflows
    
def create_graph(model, attr:str = 'mkt_share', timestamp:str = None):
    if timestamp is None:
        timestamp = np.random.choice(data_mktshare.index, 1)[0]

    edgelist = pd.DataFrame({'source':CPs}).merge(pd.DataFrame({'target':CPs}), how = 'cross')
    G = nx.from_pandas_edgelist(edgelist, source='source', target='target')

    assert attr in ['mkt_share', 'mkt_share_delta', 'inflow', 'outflow']
    if attr == 'mkt_share':
        nx.set_node_attributes(G, data_mktshare[data_mktshare.index == timestamp].reset_index(drop=True).T.to_dict()[0], attr)
    elif attr == 'mkt_share_delta': 
        nx.set_node_attributes(G, data_delta_mktshare[data_delta_mktshare.index == timestamp].reset_index(drop=True).T.to_dict()[0], attr)
    elif attr in ['inflow', 'outflow']:
        nx.set_node_attributes(G, data_delta_mktshare[data_delta_mktshare.index == timestamp].reset_index(drop=True).T.to_dict()[0], 'mkt_share_delta')
        if attr == 'inflow':
            vals = prepare_inflows(model)
        elif attr == 'outflow':
            vals = prepare_outflows(model)

        vals = vals[vals.timestamp == timestamp].drop('timestamp', axis = 1)
        kvals = list(vals.loc[:,['source','target']].itertuples(index=False, name=None))
        vvals = vals[attr].values.tolist()
        nx.set_edge_attributes(G, dict(zip(kvals, vvals)), attr)            

    COLOR_SCHEME = "RdBu"
    colors = [G.nodes[node]['mkt_share_delta' if attr in ['inflow', 'outflow'] else attr] for node in list(G.nodes())]
    if attr in ['inflow', 'outflow']:
        norm = plt.Normalize()
        edge_attr = nx.get_edge_attributes(G, attr)
        edge_colors = plt.cm.RdBu(norm(list(edge_attr.values())))

        nx.draw_networkx(G, node_color=colors, cmap=COLOR_SCHEME, edge_color = edge_colors)#, edges = list(edge_attr.keys()))

    else:
        nx.draw_networkx(G, node_color=colors, cmap=COLOR_SCHEME)#, vmin = -1, vmax = 1)

    ax=plt.gca()
    PCM=ax.get_children()[0]  
    plt.colorbar(PCM, ax=ax)
    plt.show()
    return vals
In [16]:
create_graph(model, 'inflow')
Out[16]:
inflow target source
183 0.044975 CP_1 CP_7
195 0.102212 CP_2 CP_0
279 0.052960 CP_2 CP_7
423 0.033938 CP_4 CP_3
447 0.042164 CP_4 CP_5
471 0.014331 CP_4 CP_7
639 0.034018 CP_6 CP_5

Conclusiones

En este artículo hemos explorado cómo asignar los flujos de cuota de mercado trasvasados entre distintos competidores, mediante un modelo de programación lineal escrito con pyomo y que hemos resuelto con solvers open source, y cuyos resultados hemos visualizado en forma de grafo recurriendo a la librería networkx.

En lo referente al modelado, además de variables, conjuntos y restricciones habituales en la literatura de la programación lineal, hemos usado restricciones de tipo big M así como una linealización de las restricciones de mínimo y máximo para los flujos de entrada y de salida. Como siguientes pasos, podríamos destacar que una vez determinados dichos flujos de trasvase de cuota de mercado, podríamos expandir nuestro modelo, para repercutirlo a parámetros como el diferencial económico entre los distintos competidores, acciones de mercado o percepción en plataformas como RRSS.

Información de sesión

In [17]:
import session_info
In [18]:
session_info.show(html=False)
-----
matplotlib          3.6.2
networkx            3.0
numpy               1.23.4
pandas              1.4.4
pyomo               6.4.4
session_info        1.0.0
tqdm                4.64.1
-----
IPython             8.6.0
jupyter_client      7.4.7
jupyter_core        5.2.0
jupyterlab          3.3.4
notebook            6.5.2
-----
Python 3.8.15 (default, Nov  4 2022, 15:16:59) [MSC v.1916 64 bit (AMD64)]
Windows-10-10.0.19044-SP0
-----
Session information updated at 2023-01-31 20:33

Bibliografía


Proyecto COIN-OR web
Linealización de restricciones de mínimos web
Linealización de restricciones de máximos web

¿Cómo citar este documento?

Uso de la programación lineal para la resolución de un modelo de flujos de cuota de mercado by Francisco Espiga, available under a CC BY-NC-SA 4.0 at https://www.cienciadedatos.net/documentos/py47-uso-programacion-lineal-flujos-marketshare.html DOI


¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊


Creative Commons Licence
This work by Francisco Espiga is licensed under a Creative Commons Attribution 4.0 International License.