数学建模_饮食计划

摘要：

民以食为天，合理的饮食是身体健康的基础。科学的控制摄入食物的比例可以健康的减肥。实际的饮食计划中既要考虑较低的热量摄入，还要考虑较高的满足感和饱腹感，并且营养要均衡。本文采用多目标加权、分优先级的方法将多目标优化问题转化为多个线性规划问题并利用scipy求解。最终得出一份综合多个目标、符合已知约束条件的饮食表。

一：提出问题

1.1、问题描述

在不同的时候需要调整自己的饮食计划，包括增肌和减脂。
对饮食计划进行建模，要求：

（1）只考虑饮食，不考虑运动，假设基础代谢为常数。
（2）只考虑食物中的主要成分：蛋白质，脂肪，碳水化合物三种主要热量来源，考虑维生素，钠离子两种微量元素，不考虑其它微量元素。
（3）饮食计划要综合考虑：热量摄入、饱腹感（GI值）、满足感（高脂肪高糖满足感高）。
（4）考虑营养均衡。单种营养素的范围，脂肪中的饱和脂肪与不饱和脂肪，碳水化合物中膳食纤维。（不饱和脂肪和膳食纤维的热量贡献非常小）
（5）假设食物可以由营养素线性表出，没有损失

1.2、问题背景

中国营养学会2000年提出中国居民膳食能量参考摄入量指出，成年男性轻、中体力劳动者每日需要能量为2400-2700kcal；女性轻、中体力劳动者每日需要能量为2100-2300kcal。婴儿、儿童和青少年、孕妇和乳母、老年人各自的生理特点不同，能量需要也不尽相同。如果想要更精确的得出自己每天摄入的热量，可以参考基础代谢率公示（摄入量大于这个数字）：

女: BMR = 655 + ( 9.6 x 体重kg ) + ( 1.8 x 身高cm ) - ( 4.7 x 年龄years )
男: BMR = 660 + ( 13.7 x 体重kg ) + ( 5 x 身高cm ) - ( 6.8 x 年龄years )

其中蛋白质摄取量应为人体每日所需热量的10%～15%；碳水化合物摄取量应不少于人体每日所需热量的 55%；脂肪的摄取量应不超过每日所需热量的30%。此外，每天摄取的盐不应超过6克，膳食纤维每天的摄取量应不少于16克。
人类生存需要能量，并从食物中获取该能量。食物中的卡路里含量是该食品产生多少潜在能量的量度标准。食物一般由碳水化合物、蛋白质、脂肪这三种物质组成。因此只要知道食物中这三种物质的含量，就可以知道食物含多少卡路里或多少能量。

1、碳水化合物产生热能 = 4 千卡（4大卡）/克
2、蛋白质产生热量 = 4 千卡（4大卡）/克
3、脂肪产生热量 = 9 千卡（9大卡）/克

维生素C支持人体的保卫系统和缔结组织。它在运动时会随汗液排出体外，必须及时补充。维生素E保护人身不受外来恶劣环境的侵害，尤其是过度劳累和紧张，环境破坏及训练等对人体的侵害。同时促进肌肉对氧的吸收。要满足人体对维生素的需要，只有通过食用富含维生素的新鲜食物，如水果、蔬菜和全麦制品。对于那些平时维生素摄入不足，膳食不平衡的人，可服用维生素补剂。成人一天需要摄入的各种维生素的量如下：维生素A原或前维生素A 2毫克；维生素B1 1.4毫克；维生素B2 1.6毫克；维生素B6 2毫克；维生素B12 1微克；维生素C 60毫克；维生素E 10毫克；尼克酸 18毫克；叶酸 200微克；泛酸 6毫克。

人体内钠在一般情况下不易缺乏、但在某些情况下，如禁食、少食，膳食钠限制过严而摄入非常低时，或在高温、重体力劳动、过量出汗、肠胃疾病、反复呕吐、腹泻使钠过量排出而丢失时，或某些疾病，如艾迪生病引起肾不能有效保留钠时，胃肠外营养缺钠或低钠时，利尿剂的使用而抑制肾小管重吸收钠时均可引起钠缺乏。钠的缺乏在早期症状不明显，倦怠、淡漠、无神、甚至起立时昏倒。失钠达0.5g/kg体重以上时，可出现恶心、呕吐、血压下降、痛性吉尔痉挛，尿中无氯化物检出。18岁以后建议的钠摄入量为2000毫克。

综合上述材料得出以下表格：

蛋白质脂肪(饱和脂肪) 碳水化合物(除膳食纤维) 维生素(VC) 钠离子
热量占比 %10-%15 <30% >55% 0 0
推荐摄入量热量4 热量9 热量*4 60mg 2000mg
饱腹感 2 1 2 0 0
满足感 1 2 2 0 0

中国居民平衡膳食宝塔是根据中国居民膳食指南结合中国居民的膳食结构特点设计的。它把平衡膳食的原则转化成各类食物的重量，并以直观的宝塔形式表现出来，便于群众理解和在日常生活中实行。

平衡膳食宝塔提出了一个营养上比较理想的膳食模式。它所建议的食物量，特别是奶类和豆类食物的量可能与大多数人当前的实际膳食还有一定距离，对某些贫困地区来讲可能距离还很远，但为了改善中国居民的膳食营养状况，这是不可缺的。应把它看作是一个奋斗目标，努力争取，逐步达到。
在这里插入图片描述

二：做出假设

题目中已有的假设：

只考虑饮食，不考虑运动，假设基础代谢为常数。
只考虑食物中的主要成分：蛋白质，脂肪，碳水化合物三种主要热量来源。
考虑维生素，钠离子两种微量元素，不考虑其它微量元素。
假设食物可以由营养素线性表出，没有损失

关于代谢

人一天的基础代谢可由BMR公式得出
基础代谢影响因素只考虑年龄、性别、身高和体重

关于营养素：

不饱和脂肪产生的热量忽略不计
膳食纤维产生的热量忽略不计
维生素只考虑维生素C一种，其他微量的忽略不计。
热量、饱腹感、满足感均只来源于蛋白质、脂肪、碳水化合物三类物质
摄入每克蛋白质、脂肪、碳水化合物得到的热量为4、9、4 kj
摄入每克蛋白质、脂肪、碳水化合物得到的饱腹感比例为2：1：2
摄入每克蛋白质、脂肪、碳水化合物得到的满足感比例为1：2：2

关于食物：

食物只考虑列出的11种
不考虑同类食物之间的差异
食物中含有的营养素以表格中给出的为准

三：建模过程

3.1、总体思路

本题要在已知的约束条件下，摄入最少的热量，其次是要获得较高的满足感和饱腹感。

分析题目可知，本题是多目标优化问题，涉及到三个目标：摄入最少的热量与获得较高的满足感和饱腹感，并且要满足已知的约束条件。

这里对三个目标的处理方法是：给三个目标分优先级依次处理，上一次处理的结果施加一定的权值作为下一次的约束条件。这样就把此多目标优化问题转化成了多个线性规划问题，然后求解。

但是由于三者的数值不具有很强的可比性，对三者赋予的权值难以确定，所以这里尝试各种比例并做出变化图，最终确定一个符合要求的权值比例并得出结果。

3.2、数据处理

首先根据人物信息(性别、年龄、身高、体重)计算出BMR数值，该值作为静态的基础代谢值。

然后打开数据文件diet_plan.csv，从中读取所有的食物中各种营养素的占比(百分比)，以小数的形式给出，读取的结果存到字典里面，字典名是食物名，字典的键是营养素名，键对应的值是该营养素在该食物中的占比(例如：milk[protain]=k,表示每克牛奶中有k克的蛋白质)。这里是动态的读取，即在运行程序前没有指定食物的名称也就是没有指定字典的名称，实现的方法是利用locals获取当前局部变量然后修改(用eval函数)。

现在就获取到了所有食物的列表和一系列储存着食物中各营养素比例的字典。
在这里插入图片描述

3.3、建模

决策变量：

这里的决策变量是各种食物每天的摄入量

目标函数：

这里的目标函数有三个：最小化热量、最大化GI、最大化满足感

热量：
热量来源有蛋白质、脂肪、碳水三个部分。
计算公式是每克蛋白质、脂肪、碳水产能4、9、4kj
饱腹感(GI)与满足感：
饱腹感与满足感同样来源于蛋白质、脂肪、碳水三个部分。
查阅资料可知蛋白质、碳水饱腹感强于脂肪，所以自定义三者单位质量饱腹感比例为2：1：1。
而脂肪与碳水的满足感要高于蛋白质，所以自定义三者单位质量满足感比例为1：2：2。
这里直接用摄入食物的质量乘以比例，得到的数值没有具体含义，数值的变化趋势才有意义，后面作图时对其做了一定的放大方便将三者放在一起观察。

1.# 三个目标函数
2.# min all_heat
3.c_heat = np.sum([vec_protain * 4, vec_fat * 9, vec_sugar * 4], axis=0)
4.# max gi
5.c_gi = np.sum([vec_protain * 2, vec_fat * 1, vec_sugar * 2], axis=0)
6.c_gi *= -1.5  # 符号是由于目标是最大化，数值是用来放大方便再图上观察(数值无意义、变化与比例才有意义)
7.# max satisfaction
8.c_sat = np.sum([vec_protain * 1, vec_fat * 2, vec_sugar * 2], axis=0)
9.c_sat *= -1.8  # 同上

约束条件：

1.摄入的热量要高于BMR
2.蛋白质摄取量应为人体每日所需热量的10%～15%；
3.碳水化合物摄取量应不少于人体每日所需热量的 55%；
4.脂肪的摄取量应不超过每日所需热量的30%。
5.维生素摄入60mg左右
6.钠离子摄入2000mg左右

# 约束条件，A[i]*x[i] < b[i]
A = [vec_protain * 4, -4 * vec_protain, vec_fat * 9, -9 * vec_fat, -4 * vec_sugar, 4 * vec_sugar, vec_vita, -1 * vec_vita, vec_na, -1 * vec_na, -1 * c_heat]
b = [0.15 * BMR, -0.1 * BMR, 0.3 * BMR, 0.1 * BMR, -0.55 * BMR, 0.8 * BMR, 8.8, -7.2, 7, -5, -1 * BMR]

3.4、Python求解

这里用SciPy.linprog求解,该模块标准形式如下：
在这里插入图片描述
这里对三个目标的处理方法是：给三个目标分优先级依次处理，上一次处理的结果施加一定的权值作为下一次的约束条件。这样就把此多目标优化问题转化成了多个线性规划问题，然后求解。但是由于三者的数值不具有很强的可比性，对三者赋予的权值难以确定，所以这里尝试各种比例并做出变化图，最终确定一个符合要求的权值比例并得出结果。

低热量、高GI、高满足感三者占比为1：i：j。其中ij均以step为步长从beg增长到end。首先按照初始约束计算优先级最低的满足感，然后依据其算得的结果拓展约束条件(后续得到的满足感不低于此满足感i)，然后再计算优先级中等的饱腹感，再依据其算得的结果拓展约束条件(后续得到的饱腹感不低于此饱腹感j)，最后根据此时的约束条件计算优先级最高的热量摄入。

通过这样的方法就把几种目标函数以不同的权值组合在一起，结果以3D图片的形式展示出来(见下图，xy分别为占比，z轴对于三个目标函数的度量)，并把计算过程中的数据储存在data.txt中，运算过程中如果没有正确的规划出方案，那么这一次的方案就取上一次的结果，并且写入错误日志log.txt中。

在这里插入图片描述

1.for i in np.arange(beg, end, step):
2.    for j in np.arange(beg, end, step):
3.        A = [vec_protain * 4, -4 * vec_protain, vec_fat * 9, -9 * vec_fat, -4 * vec_sugar, 4 * vec_sugar, vec_vita, -1 * vec_vita, vec_na, -1 * vec_na, -1 * c_heat]
4.        b = [0.15 * BMR, -0.1 * BMR, 0.3 * BMR, 0.1 * BMR, -0.55 * BMR, 0.8 * BMR, 8.8, -7.2, 7, -5, -1 * BMR]
5.
6.        c = c_sat  # 首先是优先级最低的: 满足感
7.        res = linprog(c, A, b, None, None, None)  # 进行计算
8.        if res.success:
9.            A.append(c)  # 加上约束条件
10.            b.append(res.fun * i)  # 加上权值i，然后添加到约束条件
11.
12.            c = c_gi  # 然后是优先级其次的饱腹感
13.            res = linprog(c, A, b, None, None, None)  # 进行计算
14.            if res.success:
15.                A.append(c)  # 加上约束条件
16.                b.append(res.fun * j)  # 加上权值j，然后添加到约束条件
17.
18.                c = c_heat  # 最后是优先级最高的满足感
19.                res = linprog(c, A, b, None, None, None)  # 进行计算

分析上述结果，最终三者的比例取值为100：95：95，即低热量占比略高于其他两个目标。
然后按照此比例进行计算(最终代码内将这部分放到了前面避免后续计算时间过长影响使用，对结论无影响)

for c in [c_sat, c_gi, c_heat]:  # 优先级排序：低热量>饱腹感>满足感
    res = linprog(c, A, b, None, None, None)  # 进行计算

    A.append(c)  # 加上约束条件
    b.append(res.fun * 0.95)  # 损失不超过5%，三者占比95：95：100

四：结果分析

三次添加约束条件计算的中间结果如下：
在这里插入图片描述
可见第一次以最大化满足感获得的结果分别是2904，2108，1771.这里的满足感是三次最高的。然后将满足感作为条件加入已有的约束，并以最大化饱腹感为目标获得的结果是2904，2108，1771，再将其加入约束，计算最小化摄入量，得到的结果为2641，2003，1727。可见摄入量较前两次有明显的减少，但代价是饱腹感和满足感的减少，且减少的幅度都在约束范围内。
所以最终的结果是最小摄入量2641kj，最大饱腹感1727，最大满足感2003。
达到此效果饮食的各成分含量如下：

在这里插入图片描述
从上图可见，欲达到上述效果，每天饮食成分为：
食用油4g，盐5g，奶类 45g，大豆及坚果 7g，鱼 80g，畜禽肉 8g，蛋类 7g，蔬菜类 157g，水果类 220g，全谷杂豆类 220g，薯类 311g

五：问题与不足

数据的准确性：

很多数据在网上查不到权威的资料，导致实验数据量小且数据不够全面，会导致得出的结论与实际情况有所出入。

实验方法：

实验将多目标通过优先级和加权的方法转换成线性规划问题，虽然一定程度上简化并解决了问题，但是也忽略了多个之间的关系。这样得出的结果有较大的局限性。
实验中选取权值的方法是利用实验中计算出的数据和生产的图片人工分析所有权值的组合，找到较合适的一组作为最后产生实验结论的权值。这样的方法不够准确，可能由于人为判断的失误导致实验的误差。

附录一：数据

,protein,fat,sugar,vitamin,na
oil,0,1,0,0,0
salt,0,0,0,0,1
milk,0.03,0.037,0.048,0,0.006
soybean_nut,0.35,0.3,0.1,0,0
fish,0.2,0.05,0.01,0,0
meat,0.18,0.25,0.1,0,0
egg,0.16,0.3,0.03,0.01,0
vegetable,0.06,0.01,0.1,0.01,0
fruit,0.05,0.01,0.15,0.01,0
grain,0.1,0.05,0.65,0.01,0
potato,0.05,0.05,0.8,0,0

附录二：代码

from scipy.optimize import *
import numpy as np
import matplotlib.pyplot as plt


'''****************************************数据初始化********************************************'''
sex = "M"  # M(man)
weight = 68  # 体重kg
height = 175  # 身高cm
age = 21  # 年龄years
local_names = locals()  # locals() 函数会以字典类型返回当前位置的全部局部变量,用于动态的从文件读取数据
foods = ['oil', 'salt', 'milk', 'soybean_nut', 'fish', 'meat', 'egg', 'vegetable', 'fruit', 'grain', 'potato']

if sex == "M":  # 男性BMR公式
    BMR = 660 + (13.7 * weight) + (5 * height) - (6.8 * age)
else:
    BMR = 655 + (9.6 * weight) + (1.8 * height) - (4.7 * age)

print("Begin".center(120, "*"))
print("\tsex:", sex, "\tweight:", weight, "\theight:", height, "\tage:", age, "\tBMR:", int(BMR), "\n\n")

mfile = open("diet_plan.csv", "r+")  # 打开储存信息的文件
names = mfile.readline().strip('\n').split(',')  # # .strip('\n'):去掉行末的换行符，并用逗号分开字符串

while True:
    line = mfile.readline()
    L = line.strip('\n').split(',')
    if not line:  # 读完停止
        break

    local_names[L[0]] = {names[1]: L[1], names[2]: L[2], names[3]: L[3], names[4]: L[4], names[5]: L[5]}
    print(L[0].center(20, "-"), eval(L[0]))
mfile.close()

print("****" * 30)

'''****************************************线性规划********************************************'''
# 用linprog

all_heat = np.zeros(len(foods))  # 热量
vec_protain = np.zeros(len(foods))  # 所有食物的蛋白质占比依次组成的向量
vec_fat = np.zeros(len(foods))
vec_sugar = np.zeros(len(foods))
vec_vita = np.zeros(len(foods))
vec_na = np.zeros(len(foods))

for i in range(len(foods)):  # 从字典向上述向量填值
    vec_protain[i] = float(eval(foods[i])['protein'])
    vec_fat[i] = float(eval(foods[i])['fat'])
    vec_sugar[i] = float(eval(foods[i])['sugar'])
    vec_vita[i] = float(eval(foods[i])['vitamin'])
    vec_na[i] = float(eval(foods[i])['na'])

# 三个目标函数
# min all_heat
c_heat = np.sum([vec_protain * 4, vec_fat * 9, vec_sugar * 4], axis=0)
# max gi
c_gi = np.sum([vec_protain * 2, vec_fat * 1, vec_sugar * 2], axis=0)
c_gi *= -1.5  # 符号是由于目标是最大化，数值是用来放大方便再图上观察(数值无意义、变化与比例才有意义)
# max satisfaction
c_sat = np.sum([vec_protain * 1, vec_fat * 2, vec_sugar * 2], axis=0)
c_sat *= -1.8  # 同上

# 约束条件，A[i]*x[i] < b[i]
A = [vec_protain * 4, -4 * vec_protain, vec_fat * 9, -9 * vec_fat, -4 * vec_sugar, 4 * vec_sugar, vec_vita,
     -1 * vec_vita, vec_na, -1 * vec_na, -1 * c_heat]
b = [0.15 * BMR, -0.1 * BMR, 0.3 * BMR, 0.1 * BMR, -0.55 * BMR, 0.8 * BMR, 8.8, -7.2, 7, -5, -1 * BMR]

for c in [c_sat, c_gi, c_heat]:  # 优先级排序：低热量>饱腹感>满足感
    res = linprog(c, A, b, None, None, None)  # 进行计算

    A.append(c)  # 加上约束条件
    b.append(res.fun * 0.95)  # 损失不超过5%，三者占比95：95：100
    # print(res.fun,c)

    if res.success:  # 能求出最优解的话
        for i in range(len(foods)):
            print(str(int(res.x[i])).rjust(5, "-"), end="")
        print()
        print("Heat :".rjust(22, "-"), abs(sum(res.x * c_heat)))
        print("Sat :".rjust(22, "-"), abs(sum(res.x * c_sat)))
        print("GI :".rjust(22, "-"), abs(sum(res.x * c_gi)))
    else:
        print("fail")

print("****" * 30, "\n", res, "\n", "****" * 30)
'''得出结论：当前目标下各种食物每日的摄入量'''
for i in range(len(foods)):
    print(foods[i].rjust(20, "-"), ":\t", int(res.x[i]), "g")

print("BMR :".rjust(22, "-"), BMR)
print("Heat :".rjust(22, "-"), sum(res.x * c_heat))
print("End".center(120, "*"))

# 至此计算完毕，下面对不同占比作图分析
'''****************************************作图分析********************************************'''
# xy轴分别是满足感和饱腹感的占比，z轴是热量摄入、满足感、饱腹感三者的结果
# 具体的占比是：热量：满足感：饱腹感 = 1：xx：yy

xx = []
yy = []
z_h = []
z_g = []
z_s = []

beg = 0.1
end = 0.9
step = 0.01
# xx = np.arange(beg, end, step) # 用于做曲面图
# yy = np.arange(beg, end, step)

log = open("log.txt", "w+")  # 错误日志信息
data = open("data.txt", "w+")  # 计算出的数据信息
for i in np.arange(beg, end, step):
    for j in np.arange(beg, end, step):
        xx.append(i)
        yy.append(j)

        A = [vec_protain * 4, -4 * vec_protain, vec_fat * 9, -9 * vec_fat, -4 * vec_sugar, 4 * vec_sugar, vec_vita,
             -1 * vec_vita, vec_na, -1 * vec_na, -1 * c_heat]
        b = [0.15 * BMR, -0.1 * BMR, 0.3 * BMR, 0.1 * BMR, -0.55 * BMR, 0.8 * BMR, 8.8, -7.2, 7, -5, -1 * BMR]

        c = c_sat  # 首先是优先级最低的: 满足感
        res = linprog(c, A, b, None, None, None)  # 进行计算
        if res.success:
            A.append(c)  # 加上约束条件
            b.append(res.fun * i)  # 加上权值i，然后添加到约束条件

            c = c_gi  # 然后是优先级其次的饱腹感
            res = linprog(c, A, b, None, None, None)  # 进行计算
            if res.success:
                A.append(c)  # 加上约束条件
                b.append(res.fun * j)  # 加上权值j，然后添加到约束条件

                c = c_heat  # 最后是优先级最高的满足感
                res = linprog(c, A, b, None, None, None)  # 进行计算
                if res.success:
                    data.write(str(i)+","+str(j)+","+str(abs(sum(res.x * c_heat)))+","
                               +str(abs(sum(res.x * c_sat)))+","+str(abs(sum(res.x * c_gi)))+str("\n"))
                    z_h.append(abs(sum(res.x * c_heat)))
                    z_s.append(abs(sum(res.x * c_sat)))
                    z_g.append(abs(sum(res.x * c_gi)))
                    continue
        # print(i, "__", j, "__fail")
        log.write(str(i)+","+str(j)+"\n")  # 写错误信息到log.txt
        z_h.append(z_h[-1])
        z_s.append(z_s[-1])
        z_g.append(z_g[-1])
    # 程序运行较久，故打印进度
    print("\r\t\t", int((i - beg + step) / (end - beg) * 100), "%", end="")

log.close()
data.close()
# 以曲面的方式画出
# z_h = np.array(z_h).reshape([len(xx), len(yy)]).T
# z_g = np.array(z_g).reshape([len(xx), len(yy)]).T
# z_s = np.array(z_s).reshape([len(xx), len(yy)]).T
# fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
# ax.plot_surface(np.array(xx), np.array(yy), np.array(z_h), color='red', linewidth=0, antialiased=False)
# ax.plot_surface(np.array(xx), np.array(yy), np.array(z_g), color='green', linewidth=0, antialiased=False)
# ax.plot_surface(np.array(xx), np.array(yy), np.array(z_s), color='blue', linewidth=0, antialiased=False)
'''以曲线的方式画出'''
ax = plt.axes(projection='3d')
# ax.plot3D(np.array(xx), np.array(yy), np.array(z_h), 'red')
# ax.plot3D(np.array(xx), np.array(yy), np.array(z_g), 'green')
# ax.plot3D(np.array(xx), np.array(yy), np.array(z_s), 'blue')
'''以散点图的方式画出'''
ax.scatter3D(np.array(xx), np.array(yy), np.array(z_h), c='red',s=1,label="Min Heat")
ax.scatter3D(np.array(xx), np.array(yy), np.array(z_g), c='green',s=1, label="Max GI")
ax.scatter3D(np.array(xx), np.array(yy), np.array(z_s), c='blue',s=1, label="Max Satisfaction")


ax.set_xlabel('satisfaction rate')
ax.set_ylabel('GI rate')
ax.set_zlabel('Target')

plt.legend()
plt.savefig("result.png")
plt.show()

posted @ 2022-03-05 15:59 Cheney822 阅读(376) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Cheney' Blog

Welcome To Cheney' Blog