数学建模作业

第五章5.4
学号后四位：3027

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
plt.rcParams['font.sans-serif']=['Times New Roman + SimSun + WFM Sans SC']
plt.rcParams['mathtext.fontset']='cm'
plt.rcParams['axes.unicode_minus']=False   
plt.rcParams['figure.dpi'] = 200
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
a = np.arange(100, 0, -1)
x = cp.Variable(100, nonneg=True)
obj = cp.Maximize(cp.sum(cp.sqrt(x)))
cons = [
    x[0] <= 10,
    x[0] + 2*x[1] <= 20,
    x[0] + 2*x[1] + 3*x[2] <= 30,
    x[0] + 2*x[1] + 3*x[2] + 4*x[3] <= 40,
    cp.sum(cp.multiply(a, x)) <= 1000,
]
prob = cp.Problem(obj, cons)
prob.solve(solver='ECOS')
print(f'最优解为：{x.value}'); print(f'最优值为：{prob.value}')`

5.5

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
def obj(x):
    x1, x2, x3 = x
    return (-1)*(2*x1 + 3*x1**2 + 3*x2 + x2**2 + x3)

def ineq(x):
    x1, x2, x3 = x
    return [
        10 - (x1 + 2*x1**2 + x2 + 2*x2**2 + x3),
        50 - (x1 + x1**2 + x2 + x2**2 - x3),
        40 - (2*x1 + x1**2 + 2*x2 + x3),
        x1 + 2*x2 - 1,
    ]

def eq(x):
    x1, x2, x3 = x
    return x1**2 + x3 - 2

x0 = np.random.randn(3)
cons = [
    {'type': 'ineq', 'fun': ineq},
    {'type': 'eq', 'fun': eq},
]
bd = [(0, None), (None,None), (None, None)]
ret = minimize(obj, x0, constraints=cons, bounds=bd)
print(ret)
print('-'*100)
print("最优值为：", -ret.fun)

5.7

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(a*x + b*x**2)
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)
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 = a*x1 + b*x1**2
p2 = a*x2 + b*x2**2
p3 = a*x3 + b*x3**2
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))