import casadi as ca

# Define the parameters
p_1 = 5.0
p_2 = 1.0

# Create a CasADi function for the objective
x, y, z = ca.MX.sym('x'), ca.MX.sym('y'), ca.MX.sym('z')
objective = 3*x**4 + y**3 + z**2 + 2*x*y + 4*z # 3x^4 + y^3 + z^2 + 2xy + 4z
obj_func = ca.Function('obj_func', [x, y, z], [objective])

# Create an optimization problem
opti = ca.Opti()

# Add the variables to the optimization problem
x_opti = opti.variable()
y_opti = opti.variable()
z_opti = opti.variable()

# Add the constraints to the optimization problem
opti.subject_to(8*x_opti + 4*y_opti + 2*z_opti - p_1 == 0)
opti.subject_to(p_2*x_opti + y_opti - z_opti - 1 == 0)
opti.subject_to(x_opti >= 0)
opti.subject_to(y_opti >= 0)
opti.subject_to(z_opti >= 0)

# Use the CasADi function within the optimization problem
opti.minimize(obj_func(x_opti, y_opti, z_opti))

# Choose solver
opts = {'ipopt': {'print_level': 0}}
opti.solver('ipopt', opts)

# Solve the optimization problem
sol = opti.solve()

# Get the optimal solution
optimal_x = sol.value(x_opti)
optimal_y = sol.value(y_opti)
optimal_z = sol.value(z_opti)

print("Optimal solution: x =", optimal_x, ", y =", optimal_y, ", z =", optimal_z)