%% Programación lineal.
% Copyright (c) 2022 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
%% Índice:
%%     - Introducción y planteamiento del problema
%                 - Forma canónica
%                 - Notación matricial
%%  - Definiciones
%                 - Puntos factibles y región factible
%                 - Problema factible
%                 - Función objetivo
%                 - Problemas acotados y no acotados
%                 - Punto solución
%                 - Valor mínimo
%% - Propiedades
%% - Ejemplo y ejercicio propuesto
%% - Métodos numéricos: la función linprog
%% - Ejemplos y ejercicios propuestos
%
%%
%% Introducción y planteamiento del problema
%
% Los problemas de programación lineal son un caso particular de problemas de optimización matemática con restricciones en los que todas las funciones involucradas son lineales
% y las variables independientes son no negativas. Estas propiedades hacen que los problemas de programación lineal se puedan resolver con técnicas específicas muy eficientes.
% Los problemas de programación lineal se podrían resolver aplicando las condiciones necesarias KKT. Pero, las técnicas específicas de programación lineal son mucho más eficientes.
%
%% Forma canónica:
%
%     \matrix{
% \min \quad c_1 x_1 + \cdots c_M x_M \cr
% \quad \quad \quad \quad \quad s.t.  \quad a_{1,1} x_1 + \cdots a_{1,M} x_M  \leq b_1 \quad \cr
% \quad \quad \quad \quad \vdots \cr
% \quad \quad \quad \quad \quad \quad \quad  \quad a_{K,1} x_1 + \cdots a_{K,M} x_M  \leq b_K \quad \cr
% \quad \quad \quad \quad \quad x_1, \ldots x_M \geq 0
% }
%
% donde \mathbf x =(x_1 \cdots x_M)^T son las variables independientes y \{c_m\}, \ \{a_{k,m}\}, \ \{b_k\}  son los coeficentes que definen el problema.
% Si originalmente el problema tiene alguna restricción de igualdad, esta se puede escribir como dos restricciones de desigualdad:
% g(\mathbf x) = 0  \ \Rightarrow   \ \ g(\mathbf x) \leq 0, \ -g(\mathbf x) \leq0
%
%% Notación matricial:
%
% Usando notación vectorial, el problema se puede expresar así
%             \matrix{
% \min \quad \mathbf c^T \mathbf{x} \cr
% \quad \quad \quad s.t.  \quad \mathbf a_1^T  \mathbf x  \leq b_1 \quad \cr
% \quad \quad \quad \quad \vdots \cr
% \quad \quad \quad \quad \quad \mathbf a_K^T  \mathbf x  \leq b_K \quad \cr
% \quad \quad \quad \quad \mathbf x \geq \mathbf 0
% }
% donde    \mathbf c = \left( \matrix{c_1 \cr \vdots \cr c_M} \right ), \quad
% \mathbf a_k = \left( \matrix{a_{k,1}  \cr \vdots \cr a_{k,M}} \right ), \quad k=1,\ldots,K
%
% En notación matricial las restricciones del problema se puede expresar así
%             \matrix{
% \min \quad \mathbf c^T \mathbf{x} \cr
% \quad \quad \quad s.t.  \quad\mathbf A \mathbf x  \leq \mathbf b \quad \cr
% \quad \quad \quad \quad \mathbf x \geq \mathbf 0
% }
%
% donde    \mathbf A = \left( \matrix{\mathbf a_1^T  \cr \vdots \cr \mathbf a_K^T} \right ) =
% \left( \matrix{
% a_{1,1} & \cdots & a_{1,M} \cr
% \vdots & \ddots & \vdots\cr
% a_{K,1} & \cdots & a_{K,M}
% } \right ), \quad \quad
% \mathbf b = \left( \matrix{b_1 \cr \vdots \cr b_K} \right ), \quad
%
% Por tanto, un problema de programación lineal está definido por la matriz \mathbf A y por los vectores \mathbf c y \mathbf b.
%
%% Definiciones
%
% Puntos factibles:  Son los puntos que cumplen las restricciones del problema.
% Región factible:    Es la región de \mathbb{R}^M que contiene los puntos factibles.
%                 - En el caso M=2, la región factible es un polígono convexo.
%                 - En el caso M=3, la región factible es un poliedro convexo.
%                 - En el caso M>3, la región factible es un hiperpoliedro convexo.
% Problema factible:  Es un problema en el que la región factible contiene algún punto.
% Función objetivo:  Es la función a minimizar,  \mathbf c^T \mathbf x.
% Problema no acotado:  Es un problema en el que la función objetivo puede tomar valores infinitamente pequeños (- \infty) en los puntos factibles.
% Problema acotado:  Es un problema en el que la función objetivo no toma valores infinitamente pequeños (- \infty) en los puntos factibles.
% Punto solución:  Un punto solución de un problema acotado y factible es un punto factible en que la función objetivo toma el mínimo valor de entre todos los puntos factibles.
% Valor mínimo:  El valor mínimo de un problema acotado y factible es el valor de la función objetivo en los puntos solución.
%
%% Propiedades:
%
% Si la región factible es el conjunto vacío, el problema no es factible y no tiene solución.
% Si un problema factible y acotado tiene una única solución, esta es un vértice de la región factible.
% Un problema factible y acotado puede tener infinitas soluciones. En este caso las soluciones son los puntos de una arista de la región factible que conecta dos vértices adyacentes. Estos vértices también son solución.
% Las soluciones de un problema factible y acotado saturan alguna restricción. A las restricciones saturadas se las llama restricciones activas.
% En problemas factibles y acotados no existen mínimos locales que no sean mínimos globales.
%
%
%% Ejemplo 1
%
% Un ganadero de ganado vacuno puede elegir entre dos tipos de piensos compuestos, P1 y P2, que se venden en el mercado con los siguientes precios:
% P1:  0.5  euros/Kg.
% P2:  0.2 euros/Kg
% Los piensos tienen 2 tipos de nutrientes: N1 y  N2.
% Un kilo de P1 tiene 6 unidades de N1 y 14 unidades de N2
% Un kilo de P2 tiene 4 unidades de N1, 4 unidades de N2
% El ganadero quiere suministrar diariamente a cada vaca al menos
% 60 unidades de N1
% 84 unidades de N2
% Determine la cantidad de pienso de cada tipo, por vaca y dia, con el objetivo de minimizar el coste total cumpliendo los requerimientos de nutrientes.
% Detetermine también el coste total resultante por vaca y día.
%
%% Solución:
%
% Variables independientes:
%             x_1: cantidad (en Kg) de P1 por vaca y día
%             x_2: cantidad (en Kg) de P2 por vaca y día
% Función de coste:
%              \mathbf c^T \mathbf x =0.5 x_1 + 0.2 x_2
%             Por tanto,
%             \mathbf c = \left( \matrix{0.5 \cr 0.2} \right )
% Restricciones:
%             6 x_1+ 4 x_2 \geq 60 \\
% 14 x_1+ 4 x_2 \geq 84 \\
% x_1, x_2 \geq 0            \mathbf a_1^T \mathbf x = -6 x_1- 4 x_2 \leq -60 = b_1\\
% \mathbf a_2^T \mathbf x = -14 x_1- 4 x_2 \leq -84 = b_2
%             Por tanto,
%             \mathbf A = \left( \matrix{
% -6 & -4 \cr
% -14 & -4
% } \right ) ,  \quad \quad
% \mathbf b = \left( \matrix{-60 \cr -84 } \right )
% Dibuja la función de coste y la región factible:

c = [0.5 0.2]';
b = [-60 -84]';
A = [-6 -4; -14 -4];

funcion = @(x1,x2) c(1)*x1 + c(2)*x2;
restriccion1 = @(x1,x2) A(1,1)*x1 + A(1,2)*x2 - b(1);
restriccion2 = @(x1,x2) A(2,1)*x1 + A(2,2)*x2 - b(2);
restriccion3 = @(x1,x2) x1; % x1 >=0
restriccion4 = @(x1,x2) x2; % x2 >=0
region = [-1 25 -1 25];

fcontour(funcion, region, 'LevelStep', 0.1)
hold on
fimplicit(restriccion1, region,'r')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
fimplicit(restriccion4, region,'r')
title('Funcion de coste y región factible'); colorbar; xlabel('x_1'); ylabel('x_2')
hold off

%
% Vértices factibles:
% Cada vértice está asociado a una pareja de restricciones. Hay que seleccionar los que son factibles.

D12 = [-6 -4;-14 -4]; g12=[-60 -84]';
v12=D12\g12;
disp(v12)

% En este vértice:    \mathbf a_1^T \mathbf x = b_1, \quad
% \mathbf a_2^T \mathbf x =b_2, \quad
% x_1  > 0, \quad
% x_2  > 0. Por tanto, es factible

D13 = [-6 -4;1 0]; g13=[-60 0]';
v13=D13\g13;
disp(v13)

% Este vértice:    \mathbf a_1^T \mathbf x = b_1, \quad
% \mathbf a_2^T \mathbf x > b_2, \quad
% x_1  = 0, \quad
% x_2  > 0. Por tanto, no es factible

D14 = [-6 -4;0 -1]; g14=[-60 0]';
v14=D14\g14;
disp(v14)

% Este vértice:    \mathbf a_1^T \mathbf x = b_1, \quad
% \mathbf a_2^T \mathbf x < b_2, \quad
% x_3  > 0, \quad
% x_2  = 0. Por tanto, es factible

D23 = [-14 -4; -1 0]; g23=[-84 0]';
v23=D23\g23;
disp(v23)

% Este vértice:    \mathbf a_1^T \mathbf x < b_1, \quad
% \mathbf a_2^T \mathbf x =b_2, \quad
% x_1  = 0, \quad
% x_2  > 0. Por tanto, es factible

D24 = [-14 -4; 0 -1]; g24=[-84 0]';
v24=D24\g24;
disp(v24)

% Este vértice:    \mathbf a_1^T \mathbf x > b_1, \quad
% \mathbf a_2^T \mathbf x =b_2, \quad
% x_1  > 0, \quad
% x_2  = 0. Por tanto, no es factible

D34 = [-1 0; 0 -1]; g34=[0 0]';
v34=D34\g34;
disp(v34)

% Este vértice:    \mathbf a_1^T \mathbf x > b_1, \quad
% \mathbf a_2^T \mathbf x > b_2, \quad
% x_1  = 0, \quad
% x_2  = 0. Por tanto, no es factible
%
% Valor de la función de coste en los vértice factibles:

cv12 = c'* v12;
disp(cv12)
cv14 = c'* v14;
disp(cv14)
cv23 = c'* v23;
disp(cv23)

%
% Por tanto, la solución es el punto \mathbf x^* = (3, 10.5), es decir, cada vaca debe consumir 3 Kg de P1 y 10.5 Kg de P2 por día.
% Esto supone un coste diario, por vaca, de \mathbf c^T \mathbf x^* =3.6 euros. De los cuales, 1.5 euros son en P1 y 2.1 euros son en P2.
%
% Se comprueba que las restricciones activas (saturadas) son:   \mathbf a_1^T \mathbf x^*=b_1, \ \ \mathbf a_2^T \mathbf x^*=b_2,
% y que las restricciones restantes son inactivas (no saturadas):  x_1^* > 0, \ \ x_2^* > 0.
%
% Solución aplicando las condiciones KKT
% Otra forma de resolver el problema, aunque mas laboriosa, es aplicando las condiciones KKT:
% f(\mathbf x)= \mathbf c^T \mathbf x
% h_k(\mathbf x)= - \mathbf a_k^T \mathbf x + b_k, \quad \quad k=1,\ldots, 4
% Al ser lineales las funciones, sus gradientes son vectores constantes
% \nabla f(\mathbf x)= \mathbf c
% \nabla h_k(\mathbf x)= - \mathbf a_k, \quad \quad k=1,\ldots, 4
% No existen puntos irregulares, por tanto, si hay soluciones, estas cumplen las condiciones KKT.
% Se puede comprobar que la solución, \mathbf x^* = (3, 10.5), es el único punto que cumple las condiciones KKT.
%
%
%% Ejercicio1
%
% Un ganadero de ganado vacuno puede elegir entre dos tipos de piensos compuestos, P1 y P2, que se venden en el mercado con los siguientes precios:
% P1:  0.5  euros/Kg.
% P2:  0.2 euros/Kg
% Los piensos tienen 3 tipos de nutrientes: N1, N2 y N3.
% Un kilo de P1 tiene 6 unidades de N1, 14 unidades de N2 y 6 unidaces de N3
% Un kilo de P2 tiene 4 unidades de N1, 4 unidades de N2 y 12 unidaces de N3
% El ganadero quiere suministrar diariamente a cada vaca al menos
% 60 unidades de N1
% 84 unidades de N2
% 72 unidades de N3
% minimizando el gasto en pienso.
% Plantee el ejercicio como un problema de programación lineal.
% Determine los vértices factibles del problema y el valor de la función objetivo en cada uno de ellos.
% Determine la cantidad de pienso por vaca, de cada tipo, con el objetivo de minimizar el coste total cumpliendo los requerimientos de nutrientes.
% Detetermine el coste total resultante por vaca y día.
%%
%
%% La función linprog
%
% La función linprog de Matlab resuelve problemas de programación lineal. Otros lenguajes de programación disponen de bibliotecas de funciones con herramientas equivalentes.
% En general, la llamada a la función es:  x = linprog( c, A, b, Ae, be, lb, ub ).
% Las entradas a linprog son el vector \mathbf c, la matriz \mathbf A y el vector \mathbf b, descritos anteriormente. Además, la función linprog permite incluir explicitamente restricciones lineales de igualdad, en la forma \mathbf {A_e}   \mathbf x = \mathbf {b_e}, y límites superiores e inferiores para las variables independientes mediante los vectores \mathbf {l_b} y \mathbf {u_b}. Internamente, la función reescribe el problema en la forma canónica mediante transformaciones lineales de las variables, y después aplica un algoritmo de programación lineal para resolverlo.
% A continuación se aplica la función para obtener la solución del problema del ejemplo 1:

c = [0.5 0.2]';
b = [-60 -84]';
A = [-6 -4; -14 -4];
lb = [0 0];
ub = [inf inf];
[x,cx]=linprog(c,A,b,[],[],lb,ub);
disp(x)
disp(cx)

%
%
%% Ejercicio 2
%
% Obtenga la solución del ejercicio 1 utilizando la función linprog
%
%
%% Ejercicio 3
%
% Considere el problema del ejercicio 1. Suponga que el precio de los piensos ha cambiado:
% P1 cuesta 0.3 euros/Kg
% P2 cuesta 0.4 euros/Kg
% Determine la cantidad de pienso por vaca y día, de cada tipo, con el objetivo de minimizar el coste total cumpliendo los requerimientos de nutrientes. También deteterminar el coste total resultante por vaca y día.
%
%
%% Ejercicio 4
%
% Considere el problema del ejercicio 1. Suponga ahora que en el mercado aparece un nuevo pienso, P3 con las siguientes características:.
% Precio: 0.3 euros/Kg
% Componentes nutriciones (por Kg): 4 unidades de N1, 10 unidades de N2 y 4 unidaces de N3
% Determine la cantidad de pienso por vaca y día, de cada tipo, con el objetivo de minimizar el coste total cumpliendo los requerimientos de nutrientes. También deteterminar el coste total resultante por vaca y día.
%%
%
%% Ejemplo 2
%
% Una compañía produce diariamente 2 tipos de piezas: P1 y P2.
% Las piezas se venden en el mercado a los siguientes precios:    P1: 83 euros,   P2: 108 euros
% Para fabricar las piezas se utiliza acero que se procesa con máquinaria especializada.
%                 Producir una unidad de P1 requiere 5 Kg de acero, maquinaria funcionando durante 6 minutos y 3 horas de mano de obra.
%                 Producir una unidad de P2 requiere 3 Kg de acero, máquinaria funcionando durante 5 minutos y 5 horas de mano de obra.
% El precio del acero es 2 euros/Kg. El coste de uso de la maquinaria es de 1 euro/min. El coste de mano de obra es de 10 euros/hora.
% Las recursos para producir las piezas son limitados:  Diarimente se dispone de 200 Kg de acero, la maquinaria está disponible durante 300 minutos al día y el tiempo de mano de obra disponible es de 200 horas diarias.
% ¿Cuantas unidades de cada pieza hay que fabricar al día para maximizar los beneficios de la compañía?
% Plantee el problema como un problema de programaciónn lineal
% Dibuje la función y la región factible
% Resuelva el problema gráficamente
% Compare la solución obtenida con la proporcionada por la función "linprog" de Matlab.
% Determine el beneficio diario de la empresa
% Determine la ganancia que obtiene la empresa por cada pieza P1 y P2
% Determine si sobra algún recurso disponible (acero, tiempo de maquinaria ó mano de obra) y en qué cantidad.
%
%% Solución:
%
% Variables independientes:
%             x_1: numero de unidades de P1 fabricadas cada día
%             x_2: numero de unidades de P2 fabricadas cada día
%
% Función de coste:   Es la función opuesta al beneficio
% -83 x_1 - 108 x_2 + 10 x_1 + 6 x_2 + 6 x_1 +5  x_2 + 30 x_1 + 50 x_2 =  \mathbf c^T \mathbf x,  con  \mathbf c = \left( \matrix{-37 \cr -47} \right )
% Restricciones:
%     \mathbf a_1^T \mathbf x = 5 x_1 + 3 x_2 \leq 200 \\
% \mathbf a_2^T \mathbf x = 6 x_1 + 5 x_2 \leq 300 \\
% \mathbf a_3^T \mathbf x = 3 x_1 + 5 x_2 \leq 200 \\
% x_1  \geq 0, x_2  \geq 0
%      Por tanto,
%             \mathbf A = \left( \matrix{
% 5 & 3 \cr
% 6 & 5 \cr
% 3 & 5
% } \right ) ,  \quad \quad
%
% \mathbf b = \left( \matrix{200 \cr 300 \cr 200} \right )
%
% Dibuja la función de coste y los vértices factibles:

c = [-37 -47]';
b = [200 300 200]';
A = [5 3; 6 5; 3 5];

funcion = @(x1,x2) c(1)*x1 + c(2)*x2;
restriccion1 = @(x1,x2) A(1,1)*x1 + A(1,2)*x2 - b(1);
restriccion2 = @(x1,x2) A(2,1)*x1 + A(2,2)*x2 - b(2);
restriccion3 = @(x1,x2) A(3,1)*x1 + A(3,2)*x2 - b(3);
restriccion4 = @(x1,x2) x1;
restriccion5 = @(x1,x2) x2;
region = [-5 80 -5 80];

fcontour(funcion, region, 'LevelStep', 30)
hold on
fimplicit(restriccion1, region,'r')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
fimplicit(restriccion4, region,'r')
fimplicit(restriccion5, region,'r')
title('Funcion de coste y región factible'); colorbar; xlabel('x_1'); ylabel('x_2')
hold off

% Solución:

lb = [0 0];
ub = [inf inf];
[x,cx]=linprog(c,A,b,[],[],lb,ub);
disp(x)
disp(cx)
disp(['Beneficio diario = ',num2str(-cx),' euros'])

%
% Beneficio por pieza:
% P1:    -c_1 \ x_1= 25 * 37
% P2:    -c_2 \ x_2= 25 * 47
%
% Recursos sobrantes:

A*x

%            \mathbf a_1^T \mathbf x = 5 x_1 + 3 x_2 = 200 \leq 200 \\
% \mathbf a_2^T \mathbf x = 6 x_1 + 5 x_2 = 275 < 300 \\
% \mathbf a_3^T \mathbf x = 3 x_1 + 5 x_2 = 200 \leq 200 \\
% x_1  > 0 \\
% x_2  > 0
% Por tanto, se saturan la primera y tercera restricción por lo que no sobra acero ni mano de obra.
% Pero no se satura la segunda restricción por lo que no se usa todo el tiempo de maquinaria disponible. Sobran 25 minutos de uso de maquinaria al día.
%
%% Ejercicio 5
%
% A la compañía del ejemplo 2 se le ofrece un contrato de 9000 euros para fabricar diariamente una pieza especial (P3), diferentes de P1 y P2, durante 10 dias.
% Cada pieza especial requiere 60 Kg de acero, maquinaria funcionando durante 100 minutos y 30 horas de mano de obra.
% ¿Debería la compañía aceptar el contrato?.
% A partir de qué precio de contrato resultaría rentable a la empresa la fabricación de las piezas especiales?.
%
%% Ejemplo 3
%
% Sea Z una variable aleatoria discreta que puede tomar 3 valores,  Z \in \Omega_Z = \{z_1, z_2, z_3\}, con probabilidades  p_1=P(Z=z_1), \ p_2=P(Z=z_2), \ p_3=P(Z=z_3)
% Conocido el espacio muestral \Omega_Z, determine la distribución de probabilidad que minimiza la media de Z, sujeto a que  p_1 \leq 1/2
%
% Variables independientes:
% p_1, p_2, p_3
%
% Función de coste:
% E[Z]  =  z_1 \ p_1 + z_2 \ p_2 + z_3 \ p_3,
% Por tanto, la función de coste  \mathbf c^T \mathbf p,  con  \mathbf c = \left( \matrix{
% z_1 \cr z_2 \cr z_3 } \right )
% Restricciones:
% p_i \geq 0,  \ i=1,2,3
% p_1 \leq 1/2, \ p_2 \leq 1, \ p_3 \leq 1
% p_1+ p_2 + p_3 = 1,
%
%% Solución:
%
% Caso particular:

z1 = 1; z2 = 2; z3 = 3;
c = [z1 z2 z3]';
lb = [0 0 0];
ub = [1/2 1 1];
Ae = [1 1 1]; be = 1;
[p,cp]=linprog(c,[],[],Ae,be,lb,ub);
disp(p)
disp(cp)

%
%% Ejercicio 6
%
% Una empresa tiene 2 plantas de producción (P1, P2), en distintos lugares, que fabrican un único producto. El producto se vende en 3 mercados (M1, M2 y M3), localizados en distintos lugares.
% En la planta P1 hay 90 unidades del producto y en la planta P2 hay 60 unidades.
% Se debe transportar 30 unidades del producto al mercado M1, 50 a M2 y 70 a M3.
% El coste (en euros) de transporte de una unidad de producto desde las distintas plantas a los distintos mercados son:
% P1 --> M1: 10,      P1 --> M2: 8,    P1 --> M3: 11
% P2 --> M1: 12,      P2 --> M2: 6,    P2 --> M3: 10
%
% Determinar el número de unidades del producto que se debe transportar desde cada planta hasta cada mercado minimizando el precio del transporte.
% De acuerdo a la solución obtenida, determine el número de unidades que no se transportan y en qué plantas están.
% Determine el coste resultante.
%
%% Ejercicio 7
%
% Un agricultor dispone de L carros de terreno (nota: el carro es una medida de área equivalente a 165  m^2 aproximadamente) donde plantar 3 cultivos diferentes: C_1, C_2 y C_3.
% El agricultor dispone de F Kg de fertilizante. Por otra parte, tiene que comprar pesticida que se vende a C_p euros el Kg.
% Cada carro de C_1 requiere F_1 Kg de fertilizante y P_1 Kg de pesticida, cada carro de C_2 requiere F_2 Kg de fertilizante y P_2 Kg de pesticida,
% y cada carro de C_3 requiere F_3 Kg de fertilizante y P_3 Kg de pesticida.
% La cosecha de un carro de C_1 se vende por S_1 euros, la cosecha de un carro de C_2 se vende por S_2 euros, y la cosecha de un carro de C_3 se vende por S_3 euros.
% El agricultor desea determinar el numero de carros que debe planta de cada cultivo.
%
% Plantee el problema de optimización
% Considere los siguientes datos:   L=20, \ F=5, \  C_p=3,  \ F_1=0.1,  \ F_2=0.2,  \ F_3=0.3,  \  P_1=0.1,  \ P_2=0.3,  \ P_3=0.4,  \  S_1=100,  \  S_2=150,  \  S_3=200)
% Resuelva el problema
% Determine el coste de pesticida resultante
% Determine el beneficio máximo
%
%% Ejemplo 4
%
% Este ejemplo está obtenido de la siguiente referencia:  Enrique Castillo et al, "Building and solving mathematical programming models in engineering and science". John Wiley and Sons, 2002.
%
% Una red de transporte es un sistema está formada por un conjunto de nudos conectados, a través de la cuál se transporta un "producto". Ejemplos típicos son las redes de tuberias, redes de carreteras y redes comunicaciones, por las que se transportan líquidos, vehículos, e información (bits), respectivamente. El objetivo es enviar cierta cantidad de producto (flujo) desde los nudos fuente hasta los nudos destino con el menor coste posible. Además, la red contiene casi siempre nudos intermedios donde no se produce ni se consume producto.
% Las redes de transporte se representan mediante grafos direccionales como el de la figura.
%
% En este ejemplo, la red  consta de N=4 nodos con los siguientes enlaces:  1 \rightarrow 2, \ \ 1 \rightarrow 3, \ \  1\rightarrow 4, \ \ 2\rightarrow 4, \ \ 3 \rightarrow 4.
% El nodo 1 es un nodo fuente con flujo entrante f_1=6.2. Los nudos 2 y 4 son nudos destino con flujos salientes  f_2=1.6,\ \   f_4=4.6. Por tanto, el producto que entra por el nodo 1 con flujo f_1  debe ser transportado a los nodos 2 y 4 con flujos f_2 y f_4, respectivamente. Los valores de los flujos entrantes y salientes a la red se consideran positivos.
% Las variables del problema (a optimizar) son \mathbf x = [x_{1,2}, \ x_{1,3}, \ x_{1,4},  \ x_{2,4},  \ x_{3,4}]^T, donde x_{i,j} representa el flujo entre los nodos  i \rightarrow j. Se considera que el flujo en las conexiones es siempre positivo.
% Además cada enlace entre dos nodos (conexión) tiene una capacidad máxima. En el ejemplo supondremos que las capacidades máximas son: m_{1,2}=3.5, \ \ m_{1,3}=3.8, \ \ m_{1,4}=2.5, \ \ m_{2,4}=4.0, \ \ m_{3,4}=3.5. Siempre se consideran valores positivos.
% Finalmente, hay un coste asociado al transporte a través de cada conexión. En este ejemplo supondremos que los costes son:  c_{1,2}=2.1, \ \ c_{1,3}=1.8, \ \ c_{1,4}=3.1, \ \ c_{2,4}=1.3, \ \ c_{3,4}=1.5. Siempre se consideran valores positivos.
%
%% Solución:
%
% A continuación se plantea el problema como un problema de programación lineal y se resuelve.
%
% Variables independientes:   \mathbf x = [x_{1,2}, \ x_{1,3}, \ x_{1,4},  \ x_{2,4},  \ x_{3,4}]^T
%
% Función de coste:   \mathbf c^T \mathbf x
% donde   \mathbf c = [c_{1,2}, \ c_{1,3}, \ c_{1,4},  \ c_{2,4},  \ c_{3,4}]^T
%
% Restricciones asociadas a los flujos de entrada y salida en cada nodo:
% x_{1,2} + x_{1,3} + x_{1,4} = f_1
% x_{1,2} - x_{2,4} = f_2
% x_{1,3} - x_{3,4} = 0
% x_{1,4} + x_{2,4} + x_{3,4} = f_4
% \mathbf A_e = \left( \matrix{
% 1 & 1 & 1 & 0 & 0 \cr
% 1 & 0 & 0 & -1 & 0 \cr
% 0 & 1 & 0 & 0 & -1 \cr
% 0 & 0 & 1 & 1 & 1
% } \right ) ,  \quad \quad
%
% \mathbf b_e = \left( \matrix{f_1 \cr f_2 \cr 0 \cr f_4} \right )
%
% Restricciones asociadas a las capacidades máximas de las conexiones entre los nodos:
% 0 \leq x_{i,j} \leq m_{i,j},
% \mathbf l_b = [0,  0,  0,   0, 0]
% \mathbf u_b = [m_{1,2},  m_{1,3},  m_{1,4},   m_{2,4},  m_{3,4}]
% \mathbf A = [], \mathbf b = []
%
% Solución:

c = [2.1,1.8,3.1,1.3,1.5]';
be = [6.2,1.6,0,4.6]';
Ae = [1,1,1,0,0; 1,0,0,-1,0; 0,1,0,0,-1; 0,0,1,1,1];
b=[];
A=[];
lb=[0,0,0,0,0];
ub=[3.5,3.8,2.5,4.0,3.5];

[x,cx]=linprog(c,A,b,Ae,be,lb,ub);
disp(x)
disp(cx)

%
% Interpretación de la solución:
% El flujo entrante por el nodo 1 (fuente) se distribuye entre los nodos 2, 3 y 4. Al nodo 2 (destino) le llega el flujo desde el nodo 1 y no suministra flujo al nodo 4. Al nodo 4 (destino) le llega el flujo diréctamente desde el nodo 1 y a través del nodo 3 (intermedio).
%