Matlab - AYUDA Problema Optimizacion Programacion Lineal

AYUDA PARA RESOLVER ESTE PROBLEMA con MATLAB (octave)

Considera el siguiente problema de optimización:

- Una empresa tiene n=2 plantas de producción, sus mercados están en m=2 ciudades destino y fabrica 1 tipo de producto.
- Cada planta tiene una capacidad máxima de producción, O(i), siendo i la planta.
- Se conocen las demandas máximas de cada producto en cada mercado, D(j), siendo j el mercado. La producción ha de ser tal que no se superen estas demandas máximas.
- El coste de transporte por unidad de producto entre la planta i y la ciudad j es C(i,j). El precio de venta del producto es el doble del máximo coste de transporte por unidad de producto.
- El número máximo de productos que se pueden transportar en cada ruta es M=6.
- La siguiente tabla resume un caso particular de este problema.
O(1) 10 C(1,1) 10
O(2) 8 C(1,2) 9
D(1) 6 C(2,1) 8
D(2) 5 C(2,2) 7

-Se debe obtener la cantidad de productos fabricados en cada planta y la distribución de unidades transportadas de manera que el beneficio total sea máximo. (La variable de decisión puede ser un vector x de tamaño 4).
Tareas a realizar:
1. Modelar el problema de optimización de forma genérica (variable y ecuaciones) cuyo objetivo sea maximizar beneficios (ingresos-costes).
2. Definir mediante matrices el problema de optimización correspondiente. Indicar a qué familia corresponde.
3. Resuelve el problema de optimización en Matlab o con Octave.

Consideraciones prácticas para realizar el código de optimización:
a. Utilizar un paso fijo de tamaño 0.001.
b. Usar un valor inicial factible estrictamente interior. Por ejemplo para las constantes de
ejemplo, el punto inicial [0 0 0 0] es una buena opción.
c. Utilizar una constante de finalización epsilon = 0.1 y una condicion se salida en función de la norma del incremento de cada paso . Utilizar como referencia la plantilla de código:

f = [];
A = [];
b = [];
x = []

X = glpk(f,A,B,[],[],"UUUU","CC",1)

t = 0.001; %Dato dado
mu = __; %Dato Desconocido
epsilon = 0.1; %Dato dado
m = size(A,1);

Gphi = A'*(1./(b-A*x));
Hphi = A'*diag((1./(b-A*x))).^2*A;
while m/t>epsilon
while (t*f+Gphi)'*Hphi*(t*f+Gphi)>epsilon
IX = -inv(Hphi)*(t*f+Gphi);
x = x+t*IX;
Gphi = A'*(1./(b-A*x));
Hphi = A'*(diag((1./(b-A*x))).^2)*A;
end
t = t*mu;
end