[Python] 라그랑주 승수법

from sympy import symbols, Eq, solve

def optimize_production(a, b, c_A, c_B, U):
    try:
        # 정의 변수
        p_A, p_B, q_A, q_B, λ = symbols('p_A p_B q_A q_B λ')

        # 목적 함수: 이윤
        profit = p_A * q_A + p_B * q_B - c_A * q_A - c_B * q_B

        # 제약 조건: 소비자 효용
        utility_constraint = Eq(a * q_A + b * q_B, U)

        # 라그랑주 함수
        Lagrangian = profit + λ * (a * q_A + b * q_B - U)

        # 편미분
        dL_dqA = Lagrangian.diff(q_A)
        dL_dqB = Lagrangian.diff(q_B)
        dL_dλ = Lagrangian.diff(λ)

        # 방정식 설정
        eq1 = Eq(dL_dqA, 0)
        eq2 = Eq(dL_dqB, 0)
        eq3 = utility_constraint

        print("eq1:", eq1)
        print("eq2:", eq2)
        print("eq3:", eq3)

        # λ 구하기
        λ_value = solve(eq1, λ)
        if not λ_value:
            raise ValueError("Failed to solve for λ")
        λ_value = λ_value[0]
        print("λ_value:", λ_value)

        # p_B 구하기
        p_B_solution = solve(eq2.subs(λ, λ_value), p_B)
        if not p_B_solution:
            raise ValueError("Failed to solve for p_B")
        p_B_solution = p_B_solution[0]
        print("p_B_solution:", p_B_solution)

        # p_A 구하기
        p_A_solution = solve(eq1.subs(p_B, p_B_solution), p_A)
        if not p_A_solution:
            raise ValueError("Failed to solve for p_A")
        p_A_solution = p_A_solution[0]
        print("p_A_solution:", p_A_solution)

        # q_B 구하기
        q_B_solution = solve(eq3.subs(q_A, U - q_A), q_B)
        if not q_B_solution:
            raise ValueError("Failed to solve for q_B")
        q_B_solution = q_B_solution[0]
        print("q_B_solution:", q_B_solution)

        # q_A 구하기
        q_A_solution = solve(eq3.subs(q_B, q_B_solution), q_A)
        if not q_A_solution:
            raise ValueError("Failed to solve for q_A")
        q_A_solution = q_A_solution[0]
        print("q_A_solution:", q_A_solution)

        return {
            'q_A': q_A_solution,
            'q_B': q_B_solution,
            'p_A': p_A_solution,
            'p_B': p_B_solution
        }
    except Exception as e:
        print(f"An error occurred: {e}")

# 예제 값
a = 1  # 제품 A의 중요도
b = 1  # 제품 B의 중요도
c_A = 40
c_B = 45
U = 200

result = optimize_production(a, b, c_A, c_B, U)
print(result)