2023高教社数学建模国赛C题 - 蔬菜类商品的自动定价与补货决策（完整代码更新）

问题一

matlab代码：

% 读取数据，假设数据文件名为 'sales_data.csv'
data = readtable('sales_data.csv');

% 数据预处理
data.销售日期 = datetime(data.销售日期, 'InputFormat', 'yyyy-MM-dd');

% 按品类汇总销售量
category_sales = varfun(@sum, data, 'GroupingVariables', '品类', 'InputVariables', '销量千克');

% 绘制蔬菜品类销售量分布图
figure;
bar(category_sales.品类, category_sales.Fun_sales千克);
title('蔬菜品类销售量分布');
xlabel('品类');
ylabel('销量(千克)');

% 分析单品销售量分布（假设要分析品类A的单品销售量）
category_A_data = data(data.品类 == '品类A', :);
single_item_sales = varfun(@sum, category_A_data, 'GroupingVariables', '单品编码', 'InputVariables', '销量千克');

% 绘制品类A单品销售量分布图
figure;
bar(single_item_sales.单品编码, single_item_sales.Fun_sales千克);
title('品类A单品销售量分布');
xlabel('单品编码');
ylabel('销量(千克)');

python代码：

import pandas as pd
import matplotlib.pyplot as plt


# 假设数据文件名为 'sales_data.csv'
data = pd.read_csv('sales_data.csv')


# 删除不需要的列
data = data.drop(columns=['扫码销售时间'])

# 转换销售日期列为日期时间格式
data['销售日期'] = pd.to_datetime(data['销售日期'])

# 可以按需要进行其他数据清洗和处理

category_sales = data.groupby('品类')['销量(千克)'].sum()


category_sales.plot(kind='bar', figsize=(10, 6))
plt.title('蔬菜品类销售量分布')
plt.xlabel('品类')
plt.ylabel('销量(千克)')
plt.show()


# 假设你想分析某一个品类的单品销售量分布，比如品类A
category_A_data = data[data['品类'] == '品类A']
single_item_sales = category_A_data.groupby('单品编码')['销量(千克)'].sum()

# 绘制单品销售量分布图
single_item_sales.plot(kind='bar', figsize=(12, 6))
plt.title('品类A单品销售量分布')
plt.xlabel('单品编码')
plt.ylabel('销量(千克)')
plt.show()

问题二

线性规划模型：

标准形式为：
min z=2X1+3X2+x
s.t
x1+4x2+2x3>=8
3x1+2x2>=6
x1,x2,x3>=0

上述线性规划问题Python代码

import numpy as np
from scipy.optimize import linprog

c = np.array([2, 3, 1])
A_up = np.array([[-1, -4, -2], [-3, -2, 0]])
b_up = np.array([-8, -6])

r = linprog(c, A_ub=A_up, b_ub=b_up, bounds=((0, None), (0, None), (0, None)))

print(r)

粒子群优化算法：

粒子群算法基本步骤
1 找出待优化的目标函数
2 设定种群规模大小（不会设置可直接采用下方代码的）
3 替换掉下方公式即可

%% 初始化种群  
f= @(x)x .* sin(x) .* cos(2 * x) - 2 * x .* sin(3 * x); % 函数表达式    % 求这个函数的最大值  
figure(1);ezplot(f,[0,0.01,20]);  
N = 50;                         % 初始种群个数  
d = 1;                          % 空间维数  
ger = 100;                      % 最大迭代次数       
limit = [0, 20];                % 设置位置参数限制  
vlimit = [-1, 1];               % 设置速度限制  
w = 0.8;                        % 惯性权重  
c1 = 0.5;                       % 自我学习因子  
c2 = 0.5;                       % 群体学习因子   
for i = 1:d  
    x = limit(i, 1) + (limit(i, 2) - limit(i, 1)) * rand(N, d);%初始种群的位置  
end  
v = rand(N, d);                  % 初始种群的速度  
xm = x;                          % 每个个体的历史最佳位置  
ym = zeros(1, d);                % 种群的历史最佳位置  
fxm = zeros(N, 1);               % 每个个体的历史最佳适应度  
fym = -inf;                      % 种群历史最佳适应度  
hold on  
plot(xm, f(xm), 'ro');title('初始状态图');  
figure(2)  
%% 群体更新  
iter = 1;  
record = zeros(ger, 1);          % 记录器  
while iter <= ger  
     fx = f(x) ; % 个体当前适应度     
     for i = 1:N        
        if fxm(i) < fx(i)  
            fxm(i) = fx(i);     % 更新个体历史最佳适应度  
            xm(i,:) = x(i,:);   % 更新个体历史最佳位置  
        end   
     end  
if fym < max(fxm)  
        [fym, nmax] = max(fxm);   % 更新群体历史最佳适应度  
        ym = xm(nmax, :);      % 更新群体历史最佳位置  
 end  
    v = v * w + c1 * rand * (xm - x) + c2 * rand * (repmat(ym, N, 1) - x);% 速度更新  
    % 边界速度处理  
    v(v > vlimit(2)) = vlimit(2);  
    v(v < vlimit(1)) = vlimit(1);  
    x = x + v;% 位置更新  
    % 边界位置处理  
    x(x > limit(2)) = limit(2);  
    x(x < limit(1)) = limit(1);  
    record(iter) = fym;%最大值记录  
     x0 = 0 : 0.01 : 20;  
     plot(x0, f(x0), 'b-', x, f(x), 'ro');title('状态位置变化')  
    pause(0.1)  
    iter = iter+1;  
end  
figure(3);plot(record);title('收敛过程')  
x0 = 0 : 0.01 : 20;  
figure(4);plot(x0, f(x0), 'b-', x, f(x), 'ro');title('最终状态位置')  
disp(['最大值：',num2str(fym)]);  
disp(['变量取值：',num2str(ym)]);  

问题三

粒子群模型：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def fit_fun(x):  # 适应函数
    return sum(100.0 * (x[0][1:] - x[0][:-1] ** 2.0) ** 2.0 + (1 - x[0][:-1]) ** 2.0)


class Particle:
    # 初始化
    def __init__(self, x_max, max_vel, dim):
        self.__pos = np.random.uniform(-x_max, x_max, (1, dim))  # 粒子的位置
        self.__vel = np.random.uniform(-max_vel, max_vel, (1, dim))  # 粒子的速度
        self.__bestPos = np.zeros((1, dim))  # 粒子最好的位置
        self.__fitnessValue = fit_fun(self.__pos)  # 适应度函数值

    def set_pos(self, value):
        self.__pos = value

    def get_pos(self):
        return self.__pos

    def set_best_pos(self, value):
        self.__bestPos = value

    def get_best_pos(self):
        return self.__bestPos

    def set_vel(self, value):
        self.__vel = value

    def get_vel(self):
        return self.__vel

    def set_fitness_value(self, value):
        self.__fitnessValue = value

    def get_fitness_value(self):
        return self.__fitnessValue


class PSO:
    def __init__(self, dim, size, iter_num, x_max, max_vel, tol, best_fitness_value=float('Inf'), C1=2, C2=2, W=1):
        self.C1 = C1
        self.C2 = C2
        self.W = W
        self.dim = dim  # 粒子的维度
        self.size = size  # 粒子个数
        self.iter_num = iter_num  # 迭代次数
        self.x_max = x_max
        self.max_vel = max_vel  # 粒子最大速度
        self.tol = tol  # 截至条件
        self.best_fitness_value = best_fitness_value
        self.best_position = np.zeros((1, dim))  # 种群最优位置
        self.fitness_val_list = []  # 每次迭代最优适应值

        # 对种群进行初始化
        self.Particle_list = [Particle(self.x_max, self.max_vel, self.dim) for i in range(self.size)]

    def set_bestFitnessValue(self, value):
        self.best_fitness_value = value

    def get_bestFitnessValue(self):
        return self.best_fitness_value

    def set_bestPosition(self, value):
        self.best_position = value

    def get_bestPosition(self):
        return self.best_position

    # 更新速度
    def update_vel(self, part):
        vel_value = self.W * part.get_vel() + self.C1 * np.random.rand() * (part.get_best_pos() - part.get_pos()) \
                    + self.C2 * np.random.rand() * (self.get_bestPosition() - part.get_pos())
        vel_value[vel_value > self.max_vel] = self.max_vel
        vel_value[vel_value < -self.max_vel] = -self.max_vel
        part.set_vel(vel_value)

    # 更新位置
    def update_pos(self, part):
        pos_value = part.get_pos() + part.get_vel()
        part.set_pos(pos_value)
        value = fit_fun(part.get_pos())
        if value < part.get_fitness_value():
            part.set_fitness_value(value)
            part.set_best_pos(pos_value)
        if value < self.get_bestFitnessValue():
            self.set_bestFitnessValue(value)
            self.set_bestPosition(pos_value)

    def update_ndim(self):

        for i in range(self.iter_num):
            for part in self.Particle_list:
                self.update_vel(part)  # 更新速度
                self.update_pos(part)  # 更新位置
            self.fitness_val_list.append(self.get_bestFitnessValue())  # 每次迭代完把当前的最优适应度存到列表
            print('第{}次最佳适应值为{}'.format(i, self.get_bestFitnessValue()))
            if self.get_bestFitnessValue() < self.tol:
                break

        return self.fitness_val_list, self.get_bestPosition()

if __name__ == '__main__':
    # test 香蕉函数
    pso = PSO(4, 5, 10000, 30, 60, 1e-4, C1=2, C2=2, W=1)
    fit_var_list, best_pos = pso.update_ndim()
    print("最优位置:" + str(best_pos))
    print("最优解:" + str(fit_var_list[-1]))
    plt.plot(range(len(fit_var_list)), fit_var_list, alpha=0.5)

模拟退火模型：

#本功能实现最小值的求解#

from matplotlib import pyplot as plt
import numpy as np
import random
import math
plt.ion()#这里需要把matplotlib改为交互状态

#初始值设定
hi=3
lo=-3
alf=0.95
T=100

#目标函数
def f(x):
    return 11*np.sin(x)+7*np.cos(5*x)##注意这里要是np.sin

#可视化函数（开始清楚一次然后重复的画）
def visual(x):
    plt.cla()
    plt.axis([lo-1,hi+1,-20,20])
    m=np.arange(lo,hi,0.0001)
    plt.plot(m,f(m))
    plt.plot(x,f(x),marker='o',color='black',markersize='4')
    plt.title('temperature={}'.format(T))
    plt.pause(0.1)#如果不停啥都看不见

#随机产生初始值
def init():
    return random.uniform(lo,hi)

#新解的随机产生
def new(x):
    x1=x+T*random.uniform(-1,1)
    if (x1<=hi)&(x1>=lo):
        return x1
    elif x1<lo:
        rand=random.uniform(-1,1)
        return rand*lo+(1-rand)*x
    else:
        rand=random.uniform(-1,1)
        return rand*hi+(1-rand)*x

#p函数
def p(x,x1):
    return math.exp(-abs(f(x)-f(x1))/T)

def main():
    global x
    global T
    x=init()
    while T>0.0001:
        visual(x)
        for i in range(500):
            x1=new(x)
            if f(x1)<=f(x):
                x=x1
            else:
                if random.random()<=p(x,x1):
                    x=x1
                else:
                    continue
        T=T*alf
    print('最小值为：{}'.format(f(x)))

main()

遗传算法模型：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

DNA_SIZE = 24
POP_SIZE = 200
CROSSOVER_RATE = 0.8
MUTATION_RATE = 0.005
N_GENERATIONS = 50
X_BOUND = [-3, 3]
Y_BOUND = [-3, 3]


def F(x, y):
	return 3*(1-x)**2*np.exp(-(x**2)-(y+1)**2)- 10*(x/5 - x**3 - y**5)*np.exp(-x**2-y**2)- 1/3**np.exp(-(x+1)**2 - y**2)

def plot_3d(ax):

	X = np.linspace(*X_BOUND, 100)
	Y = np.linspace(*Y_BOUND, 100)
	X,Y = np.meshgrid(X, Y)
	Z = F(X, Y)
	ax.plot_surface(X,Y,Z,rstride=1,cstride=1,cmap=cm.coolwarm)
	ax.set_zlim(-10,10)
	ax.set_xlabel('x')
	ax.set_ylabel('y')
	ax.set_zlabel('z')
	plt.pause(3)
	plt.show()


def get_fitness(pop): 
    x,y = translateDNA(pop)
	pred = F(x, y)
	return (pred - np.min(pred)) + 1e-3 #减去最小的适应度是为了防止适应度出现负数，通过这一步fitness的范围为[0, np.max(pred)-np.min(pred)],最后在加上一个很小的数防止出现为0的适应度


def translateDNA(pop): #pop表示种群矩阵，一行表示一个二进制编码表示的DNA，矩阵的行数为种群数目
	x_pop = pop[:,1::2]#奇数列表示X
	y_pop = pop[:,::2] #偶数列表示y
	
	#pop:(POP_SIZE,DNA_SIZE)*(DNA_SIZE,1) --> (POP_SIZE,1)
	x = x_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(X_BOUND[1]-X_BOUND[0])+X_BOUND[0]
	y = y_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(Y_BOUND[1]-Y_BOUND[0])+Y_BOUND[0]
	return x,y

def crossover_and_mutation(pop, CROSSOVER_RATE = 0.8):
	new_pop = []
	for father in pop:		#遍历种群中的每一个个体，将该个体作为父亲
		child = father		#孩子先得到父亲的全部基因（这里我把一串二进制串的那些0，1称为基因）
		if np.random.rand() < CROSSOVER_RATE:			#产生子代时不是必然发生交叉，而是以一定的概率发生交叉
			mother = pop[np.random.randint(POP_SIZE)]	#再种群中选择另一个个体，并将该个体作为母亲
			cross_points = np.random.randint(low=0, high=DNA_SIZE*2)	#随机产生交叉的点
			child[cross_points:] = mother[cross_points:]		#孩子得到位于交叉点后的母亲的基因
		mutation(child)	#每个后代有一定的机率发生变异
		new_pop.append(child)

	return new_pop

def mutation(child, MUTATION_RATE=0.003):
	if np.random.rand() < MUTATION_RATE: 				#以MUTATION_RATE的概率进行变异
		mutate_point = np.random.randint(0, DNA_SIZE)	#随机产生一个实数，代表要变异基因的位置
		child[mutate_point] = child[mutate_point]^1 	#将变异点的二进制为反转

def select(pop, fitness):    # nature selection wrt pop's fitness
    idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,
                           p=(fitness)/(fitness.sum()) )
    return pop[idx]

def print_info(pop):
	fitness = get_fitness(pop)
	max_fitness_index = np.argmax(fitness)
	print("max_fitness:", fitness[max_fitness_index])
	x,y = translateDNA(pop)
	print("最优的基因型：", pop[max_fitness_index])
	print("(x, y):", (x[max_fitness_index], y[max_fitness_index]))


if __name__ == "__main__":
	fig = plt.figure()
	ax = Axes3D(fig)	
	plt.ion()#将画图模式改为交互模式，程序遇到plt.show不会暂停，而是继续执行
	plot_3d(ax)

	pop = np.random.randint(2, size=(POP_SIZE, DNA_SIZE*2)) #matrix (POP_SIZE, DNA_SIZE)
	for _ in range(N_GENERATIONS):#迭代N代
		x,y = translateDNA(pop)
		if 'sca' in locals(): 
			sca.remove()
		sca = ax.scatter(x, y, F(x,y), c='black', marker='o');plt.show();plt.pause(0.1)
		pop = np.array(crossover_and_mutation(pop, CROSSOVER_RATE))
		#F_values = F(translateDNA(pop)[0], translateDNA(pop)[1])#x, y --> Z matrix
		fitness = get_fitness(pop)
		pop = select(pop, fitness) #选择生成新的种群
	
	print_info(pop)
	plt.ioff()
	plot_3d(ax)