import numpy as np                  # importa a biblioteca numpy
from scipy.optimize import minimize # importa a função minimize
from scipy.stats import norm        # importa a função norma

# Define a função a ser minimizada (função objetivo)
def obj(z,*args):
    beta = np.linalg.norm(z) # define beta como a norma do vetor z
    return beta              # retorna o valor de beta

# Define a restrição (g >= 0)
def cons(z,*args):
    x = mu + sigma*z # calcula o valor do vetor X a partir de Z
    g = -lsf(x)      # avalia a função de estado limite
    return g

# Define a Função de Estado Limite (g <= 0)
def lsf(x):
    b   = x[0]
    h   = x[1]
    fy  = x[2]
    q   = x[3]
    l   = x[4]
    mmax= (q*l**2)/8
    ine= (b*h**3)/12
    y=h/2
    sigma = mmax*y/ine
    g= fy-sigma #funcion de estado limite
    return g

mu    = np.array([50.0, 284.44, 250.0, 20.0, 5000.0]) # define a média das variáveis aleatórias
sigma = np.array([2.5, 10.0, 25.0, 4.0, 500.0]) # define o desvio padrão das variáveis aleatórias
z     = np.array([0.0, 0.0, 0.0, 0.0, 0.0]) # solução inicial

# Define as restrições do problema de otimização
constraints = [{'type': 'ineq', 'fun': cons}] # Inequality constraint

# Define os limites para as variáveis (e.g., -10 <= x[0] <= 10, -10 <= x[1] <= 10)
bounds = [(-10, 10), (-10, 10), (-10, 10), (-10, 10), (-10, 10)]

# Resolve o problema de otimização com a função minimize
result = minimize(obj, z, args=(mu, sigma), method='SLSQP', bounds=bounds, constraints=constraints)

z    = result.x
beta = result.fun # Nivel de confiabilidad
Pf   = norm.cdf(-beta)

print("MPP = ", z) # El punto
print("beta = ", beta)
print("Pf = ", Pf) # Probabilidad de falla
print(result.nfev) # Numero de veces que evaluo la funcion para resolver el problema