習題5.7程式碼
import cvxpy
import cvxpy as cp
import numpy as np
import pandas as pd
from scipy.optimize import minimize
import sympy as sp
sp.init_printing(use_unicode=True)
import matplotlib.pyplot as plt
x = cp.Variable(3, integer=True)
cumulative_output = cp.cumsum(x)
demand = np.array([40, 60, 80])
cumulative_demand = np.cumsum(demand)
max_output = 100
a, b, c = 50, 0.2, 4
product_fee = cp.sum(ax + bx2)
store_fee = cp.sum(c(cumulative_output - cumulative_demand)) - c(cumulative_output[-1] - cumulative_demand[-1]) # 最後一年的產量必然不大於需求
obj = cp.Minimize(product_fee + store_fee)
cons = [
cumulative_output >= cumulative_demand
]
prob = cp.Problem(obj, cons)
result = prob.solve(solver = 'GUROBI')
print(f'最優解為：{x.value}'); print(f'最優值為：{prob.value}')
x1, x2, x3 = sp.var('x1 x2 x3')
cumulative_output = [x1, x1+x2, x1+x2+x3]
a, b, c = sp.var('a b c')
p1 = ax1 + bx12
p2 = ax2 + bx22
p3 = ax3 + bx32
s1 = c(cumulative_output[0] - cumulative_demand[0])
s2 = c(cumulative_output[1] - cumulative_demand[1])
s3 = 0
p = p1 + p2 + p3
s = s1 + s2 + s3
total_cost = (p + s).simplify()
y = total_cost
dydx1 = sp.diff(y, x1)
dydx2 = sp.diff(y, x2)
dydx3 = sp.diff(y, x3)
x1min, x2min, x3min = sp.solve([dydx1, dydx2, dydx3], [x1, x2, x3]).values()
x1min, x2min, x3min
dx1minda = sp.diff(x1min, a)
dx2minda = sp.diff(x2min, a)
dx3minda = sp.diff(x3min, a)
sx1a = dx1minda * (a/x1min)
sx2a = dx1minda * (a/x2min)
sx3a = dx1minda * (a/x3min)
print('當 a=50, b=0.2, c=4 時')
print('x1對a的靈敏性：', sx1a.subs({a: 50, b: 0.2, c:4}).n(4))
print('x2對a的靈敏性：', sx2a.subs({a: 50, b: 0.2, c:4}).n(4))
print('x3對a的靈敏性：', sx3a.subs({a: 50, b: 0.2, c:4}).n(4))
dx1mindb = sp.diff(x1min, b)
dx2mindb = sp.diff(x2min, b)
dx3mindb = sp.diff(x3min, b)
sx1b = dx1mindb * (b/x1min)
sx2b = dx1mindb * (b/x2min)
sx3b = dx1mindb * (b/x3min)
print('當 a=50, b=0.2, c=4 時')
print('x1對b的靈敏性：', sx1b.subs({a: 50, b: 0.2, c:4}).n(4))
print('x2對b的靈敏性：', sx2b.subs({a: 50, b: 0.2, c:4}).n(4))
print('x3對b的靈敏性：', sx3b.subs({a: 50, b: 0.2, c:4}).n(4))
dx1mindc = sp.diff(x1min, c)
dx2mindc = sp.diff(x2min, c)
dx3mindc = sp.diff(x3min, c)
sx1c = dx1mindc * (c/x1min)
sx2c = dx1mindc * (c/x2min)
sx3c = dx1mindc * (c/x3min)
print('當 a=50, b=0.2, c=4 時')
print('x1對c的靈敏性：', sx1c.subs({a: 50, b: 0.2, c:4}).n(4))
print('x2對c的靈敏性：', sx2c.subs({a: 50, b: 0.2, c:4}).n(4))
print('x3對c的靈敏性：', sx3c.subs({a: 50, b: 0.2, c:4}).n(4))
ymin = y.subs({x1: x1min, x2: x2min, x3: x3min}).simplify()
ymin
dyminda = sp.diff(ymin, a)
dymindb = sp.diff(ymin, b)
dymindc = sp.diff(ymin, c)
sya = dyminda * (a/ymin)
syb = dymindb * (b/ymin)
syc = dymindc * (c/ymin)
print('當 a=50, b=0.2, c=4 時')
print('y對a的靈敏性：', sya.subs({a: 50, b: 0.2, c:4}).n(4))
print('y對b的靈敏性：', syb.subs({a: 50, b: 0.2, c:4}).n(4))
print('y對c的靈敏性：', syc.subs({a: 50, b: 0.2, c:4}).n(4))
