# %% Import:

import numpy as np
import matplotlib.pyplot as plt

# %% Def of functions:

def fun_f(x, y):
    f = -(4500*x + 6000*y)
    if not(y <= fun_g1(x)): f = np.inf  # Constraint
    if not(y <= fun_g2(x)): f = np.inf  # Constraint
    if not(y <= fun_g3(x)): f = np.inf  # Constraint
    return f

def fun_g1(x):
    y = -x/3 + 14/3
    return y

def fun_g2(x):
    y = -2*x/3 + 16/3
    return y

def fun_g3(x):
    y = -2*x + 12
    return y

# %% Initialization:

x_lb = 0
x_ub = 6
Nx = 1000
x_array = np.linspace(x_lb, x_ub, Nx)

y_lb = 0
y_ub = 14/3
Ny = 1000
y_array = np.linspace(y_lb, y_ub, Ny)

f_min = np.inf
x_opt = 0
y_opt = 0

# %% Grid optim:

for x in x_array:
    for y in y_array:
        f = fun_f(x, y) # Value of objective function
        if (f <= f_min):  # Improvement of solution
            f_min = f
            x_opt = x
            y_opt = y

h_max = -f_min

# %% Displaying results:

print(f'h_max = {h_max:.2f}')
print(f'x_opt = {x_opt:.2f}')
print(f'y_opt = {y_opt:.2f}')

# %% Plotting:

g1_array = fun_g1(x_array)  # Constraint fun g1(x)
g2_array = fun_g2(x_array)  # Constraint fun g2(x)
g3_array = fun_g3(x_array)  # Constraint fun g2(x)

plt.close('all')
plt.figure(1)
plt.grid()
plt.plot(x_opt, y_opt, 'ro', label='optimum')
plt.plot(x_array, g1_array, 'b', label='g1')
plt.plot(x_array, g2_array, 'g', label='g2')
plt.plot(x_array, g3_array, 'k', label='g3')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')

plt.savefig('prog_optim_lin_sofaer.pdf')
plt.show()