数学建模优秀论文分析及免费资料分享(1)——2022 E草原放牧问题 1

1. 原则与重点

本次论文与笔记会免费分享，详细视频讲解可看b站up：啊我有兔子牙，csdn下载资源不方便的也可以直接在b站视频对应评论区下载，希望与大家一起进步。

阅读一定要注意效率和时间，分析重点
每个问题的解决方法
重点方法，以及如何运用
解题的递进关系
如果是我们，到时候能怎么做 ，思考我们是专家的话，得分点是什么
可行性分析及时间分析，完成度
负责写代码的人着重看代码，搞懂输入输出

2. 方法分析

问题 1、从机理分析的角度，建立四种放牧策略如何影响锡林郭勒草原的土壤湿度和植被生物量的数学模型
问题 2、利用土壤湿度、土壤蒸发以及降水等数据，建立模型对保持目前放牧策略不变情况下对 2022 年和 2023 年共 24 个月份的 10cm、40cm、100cm 和 200cm 的土壤湿度进 行预测。
问题 3、从机理分析的角度，建立不同放牧策略（放牧方式和放牧强度）对锡林郭勒草原土壤化学性质影响的数模型。并请结合附件 14 中数据预测锡林郭勒草原监测样地 (12 个放牧小区)在不同放牧强度下 2022 年土壤同期有机碳、无机碳、全 N、土壤 C/N 比 等值。
问题 4、利用沙漠化程度指数预测模型和附件提供数据（包括自己收集的数据）确定 不同放牧强度下监测点的沙漠化程度指数值。并请尝试给出定量的土壤板结化定义，在建立合理的土壤板结化模型基础上结合问题 3，给出放牧策略模型，使得沙漠化程度指数与板结化程度最小。
问题 5、锡林郭勒草原请在给定的降水量（300mm，600mm、900 mm 和 1200mm） 四种情形下，在保持草原可持续发展情况下对锡林郭勒草原的实验草场内放牧羊的数量进 行求解，找到最大阈值。（注：这里计算结果可以不是正整数） [最优化]
问题 6、在保持附件 13 的示范牧户放牧策略不变和问题 4 中得到的放牧方案两种情况下，用图示或者动态演示方式分别预测示范区 2023 年 9 月土地状态（比如土壤肥力变化、 土壤湿度、植被覆盖等）。

针对问题一，建立不同放牧强度对土壤湿度和植被生物量影响的数学模型。
首先搜集和整合不同放牧强度下，土壤湿度和植被生物量的变化情况，对数据进行处理。根据数据和扩展阅读给出不同放牧强度对土壤湿度和植被生物量影响的微分方程模型，并结合放牧强度对不同层次土壤湿度和植被生物量影响的折线图，分别进行分析。

1. 数据进行了哪些处理
2. 为什么选择微分方程模型：
3. 进行了画图分析，但是只画了折线图

针对问题二，建立模型对 2022 年、2023 年不同深度土壤湿度进行预测。
首先对附件中的数据标准化，然后用SPSS 进行相关性分析，计算出不同土壤湿度和影响因素之间的皮尔逊相关系数，进而选取相关性显著的因素作为不同深度土壤湿度的影响因素。其次使用土壤湿度和相应的影响因素的历史数据，利用长短期记忆（LSTM）时间序列模型对 2022 年、2023 年不同深度土壤湿度进行预测。

针对问题三，建立不同放牧强度对土壤化学性质影响的数学模型，并预测在不同放牧强度下 2022 年土壤同期有机碳、无机碳、全 N、土壤 C/N 比等值。
首先将数据依据时间进行排序，并将数据定量化处理，进而标准化。其次针对每一项土壤化学性质单独设计一
个回归模型，分别调参。并画出不同放牧强度下同一营养素的箱线图，观察不同放牧强度对同一样营养素的含量影响。最后依据 2012-2020 年所提供的数据，运用 ARIMA 时间序列模型对不同放牧强度下各个放牧小区 2022 年土壤同期有机碳、无机碳、全 N、土壤 C/N 比等值进行预测。

针对问题四，建立土壤板结化模型，给出放牧策略模型。
首先对数据进行标准化，利用 R 语言 summary 函数对其做主成分分析，确定与沙漠化程度最具相关性的几个因素，并 建立表达式。其次，选取标准化后的土壤板结化相关数据作为因素做主成分分析，分别建立不同放牧强度下土壤板结化影响因素模型。最后，将第之前预测的数据带入模型，得到不同放牧强度下土壤沙漠化、板结化程度，取其最小值即为最优放牧策略。

针对问题五，在保持草原可持续发展情况下对实验草场内放牧羊的数量进行求解，找到最大阈值。
首先选取实际数据带入到问题四模型中求解参数。其次选择 0.4 和 0.3 作为沙漠化和板结化可逆的最大界限。最后使用目标规划模型，在可持续发展的前提下，利用 附件 14、15 中的数据和实际数据求解最大值，即为羊群数量的最大阈值。

针对问题六，基于时间序列模型预测 2023 年 9 月份的土壤化学性质、土壤湿度和植 被指数(NDVI)，并用图示的方式直观展示出示范区 2012-2023 年 9 月土地状态，并分析其变化趋势

3. 技术介绍

一、微分方程

微分方程是描述一个变量关于另一个变量的变化率的方程。简单来说，就是研究某个事物的变化情况是如何随着其他因素而改变的。
例如，在物理学中，牛顿第二定律 F = ma 可以写成微分方程的形式，描述了物体的加速度（变化率）与所受的力和物体质量之间的关系。在生物学中，种群增长模型也可以用微分方程来描述，比如 logistic 方程，它表示种群数量的变化率与当前种群数量以及环境容纳量等因素有关。
优点：

• 能够准确描述动态系统的变化规律，对于物理、生物等领域的复杂现象有强大的建模能力。
• 可以考虑多种因素的相互作用，通过方程中的变量和参数反映实际系统的特性。
  缺点：
• 求解复杂的微分方程通常比较困难，尤其是非线性微分方程，可能需要高级的数学方法和计算工具。
• 对于实际问题，建立准确的微分方程模型可能需要对系统有深入的了解和大量的实验数据。
  适用场景：
• 物理学中研究物体的运动、力学、电磁学等问题。
• 生物学中描述种群增长、生态系统的动态变化等。
• 工程领域中分析控制系统、电路等的动态行为。
  

二、长短期记忆（LSTM）算法

LSTM 是一种特殊的神经网络算法，主要用于处理时间序列数据，比如语音、文本、股票价格等。它能够记住长期的信息，同时也不会忘记短期的重要信息。
比如，在自然语言处理中，LSTM 可以用于机器翻译。它可以记住前面的单词和句子的信息，以便更好地翻译当前的单词和句子。在股票价格预测中，LSTM 可以根据过去一段时间的股票价格走势来预测未来的价格。
优点：

• 擅长处理时间序列数据，能够有效捕捉长期依赖关系，对于序列中的上下文信息有很好的记忆能力。
• 在自然语言处理、语音识别等任务中表现出色，可以处理变长的序列数据。
  缺点：
• 模型相对复杂，训练时间较长，需要大量的计算资源。
• 对于超大规模的数据，可能会出现过拟合的问题。
  适用场景：
• 自然语言处理任务，如机器翻译、文本生成、情感分析等。
• 语音识别和语音合成。
• 时间序列预测，如股票价格预测、气象预测等。

三、差分整合移动平均自回归模型（ARIMA）

ARIMA 是一种用于时间序列分析和预测的统计模型。它通过分析时间序列数据的过去值、当前值和随机误差项来预测未来的值。
例如，在销售预测中，可以使用 ARIMA 模型根据过去几个月的销售数据来预测未来几个月的销售量。在气象预测中，ARIMA 可以根据过去的气温、气压等数据来预测未来的天气情况。
优点：

• 理论成熟，有明确的统计性质，易于理解和解释。
• 对于平稳的时间序列数据，能够进行有效的预测。
  缺点：
• 要求数据具有平稳性，对于非平稳数据需要进行差分等预处理。
• 对于复杂的时间序列模式，可能无法准确捕捉非线性关系。
  适用场景：
• 经济领域中的销售预测、价格预测等。
• 气象数据的短期预测。
• 平稳时间序列数据的分析和预测。
  

四、主成分分析

主成分分析是一种数据降维的方法。它可以将多个相关的变量转化为少数几个不相关的综合变量，这些综合变量称为主成分。
比如，在市场调研中，可能有很多个变量来描述消费者的行为，如年龄、收入、消费习惯等。通过主成分分析，可以将这些变量转化为几个主成分，从而更方便地分析消费者的行为模式。在图像识别中，主成分分析可以用于降低图像的维度，提高识别效率。
优点：

• 可以有效地降低数据维度，减少计算量和存储空间。
• 去除数据中的冗余信息，提取主要的特征，有助于数据分析和可视化。
  缺点：
• 主成分的解释可能比较困难，尤其是在原始变量较多的情况下。
• 可能会丢失一些有用的信息，特别是对于非线性关系的描述不足。
  适用场景：
• 高维数据的降维处理，如基因表达数据、图像数据等。
• 多变量数据分析，帮助发现数据中的主要模式和趋势。
• 数据可视化，将高维数据投影到低维空间进行展示。

五、熵值法

熵值法是一种用于确定指标权重的方法。它根据指标的变异程度来确定权重，变异程度越大的指标，权重越大。
例如，在综合评价中，有多个指标来评价一个对象，如学生的学习成绩、品德、体育等。通过熵值法可以确定每个指标的权重，从而更客观地评价学生的综合表现。在企业绩效评价中，熵值法可以用于确定各个财务指标和非财务指标的权重。
优点：

• 客观地确定指标权重，避免了主观因素的影响。
• 基于数据的变异程度进行权重分配，能够反映指标的重要性。
  缺点：
• 对数据的质量要求较高，异常值可能会对结果产生较大影响。
• 只能确定指标的相对权重，不能直接反映指标的实际重要性。
  适用场景：
• 综合评价问题，如企业绩效评价、环境质量评价等。
• 多指标决策问题，确定各个指标的权重。
  案例使用：
  

六、目标规划模型

目标规划模型是一种多目标决策方法。它可以同时考虑多个目标，并在这些目标之间进行权衡和优化。
比如，在生产计划中，可能有多个目标，如生产成本最小化、产量最大化、质量最优化等。目标规划模型可以帮助企业制定最优的生产计划，在这些目标之间找到一个平衡。在资源分配问题中，目标规划模型可以根据不同的需求和约束条件，合理地分配资源。
优点：

• 能够同时考虑多个目标，并在目标之间进行权衡和优化。
• 可以处理不同类型的目标，如最大化、最小化、区间约束等。
  缺点：
• 模型的求解可能比较复杂，尤其是当目标和约束较多时。
• 对于目标之间的冲突处理可能不够灵活，需要根据具体情况进行调整。

适用场景：

• 生产计划制定，平衡成本、产量、质量等多个目标。
• 资源分配问题，满足不同部门或项目的需求。
• 多目标决策问题，在多个相互冲突的目标中寻找最优解。
  建立目标规划模型。第一步，确立目标函数即 最大的羊群数。
  第二步，写明限制条件：避免草原土壤的沙漠化或板结化到不可逆的情况 下草原的最大羊群承受能力。
  第三步，找出影响羊群数量的三个因素：植被生物量(即干重)， 放牧强度（对照 NG，轻度放牧 LGI，中度放牧 MGI，重度放牧 HGI）和湿度(土壤湿度)。
  第四步，找出三种因素与羊群数量的相关性并设表达式。第五步根据假设和数据进行模拟 预测。

七、时间序列模型
时间序列模型是专门用于分析和预测时间序列数据的模型。它可以根据过去的数据来预测未来的数据。
例如，在经济预测中，时间序列模型可以根据过去的 GDP、通货膨胀率等数据来预测未来的经济走势。在交通流量预测中，时间序列模型可以根据过去的交通流量数据来预测未来某个时间段的交通流量。
优点：

• 专门针对时间序列数据设计，能够有效地捕捉时间序列的特征和趋势。
• 可以进行短期和长期的预测，为决策提供依据。
  缺点：
• 对于复杂的时间序列模式，可能需要更复杂的模型才能准确预测。
• 对数据的质量和稳定性要求较高，异常值和数据波动可能会影响预测结果。
  适用场景：
• 经济预测、销售预测、库存管理等。
• 气象预测、交通流量预测等。
• 任何涉及时间序列数据的分析和预测问题。
  

亮点回顾（老师为什么给高分）

自评
10.1 模型优点
（1）本文通过运用统计及数学软件，例如 Python、R 语言、SPSS 等，尽可能使数据处理合理化，以得到更加准确的结果。
（2）在数据的处理方面，定性数据定量化，定量数据标准化，通过查找大量文献尽可能使数据判定标准合理化。
（3）运用多种时间序列方法对数据进行预测。长短期记忆（LSTM）时间序列模型、 差分整合移动平均自回归模型—ARIMA 模型对 2022 年之后土壤数据进行了不同类型的预 测。
（4）建立多种方程形式以减小计算误差，例如微分方程、多项式方程、一维线性方程。运用主成分分析及皮尔逊相关系数，尽可能的使用较少的因素反应更好的变量。
（5）尽可能绘制图形及表格来反映数据及模型精度，使得文章更加生动直观。

10.2 模型缺点及改进方向
（1）数据较少，许多反应问题的数据没有找到而使得模型缺少变量。后续应广泛收集数据并与相关部门合作取得更加详细得数据，从而使得建立的模型考虑得更加全面。
（2）模型方程建立多数采用一次变量，可能不能较好的拟合数据的真实情况。后续应建立真实数据拟合曲线来给自变量定阶。
（3）数据清洗不够，可能会因数据极值而使得整体数据方差较大。应先对数据进行整体观察，筛除异常值以减小异常值对模型精度的影响。
（4）由于时间原因，部分题目解答仅提供思路，对于细致的问题回答不到位。
（5）本文中涉及代码简单重复较多，没有设计出简单直接的方以直接得到所需数据。
如时间充裕，可以增强代码的一般化，使之可以解决一类问题。
客评：

1. 结构清晰内容完整充实。
   

2. 能有效拆解问题，回答问题，描述的语言简洁清晰。如：问题三一共可以分为两问、问题四可分为三个问题。

3. 论文写的很好，格式清爽，模型假设等内容很全面

4. 自己补充了建模需要的关键数据：由于附件中缺乏放牧强度对土壤物理性质影响的相关数据，我们通过广泛阅读大量文献，搜集了针对华北农牧交错地带，关于放牧强度对土壤水分影响的相关数据[3]。

5. 自己将题目中的变量转化为新的关键变量：对牧前、第一次牧后、第二次牧后、第三次牧后以及第四次牧后的各种植被干重进行求和，作为植被生物量

6. 对实验结果的分析比较深入，准确。

7. 许多处理有文献支持，数据、影响因素等，而不只是技术处理

免费资源获取

1.CSDN文章标题下方

2.csdn好像只有vip才能下载，虽然我并没设置。csdn获取麻烦的同学可以到b站up：啊我有兔子牙
工坊免费领取（自动发货需要0.01元），发货后会收到私信，资源在私信窗口打开即可。

5. 代码
   注意其中有页码被复制上去了，用的时候要删掉
   第二问

## Multivariate-SingleStep-LSTM.py 
from args import ms_args_parser 
from util import train, test 
LSTM_PATH = '../model/Multivariate-SingleStep-LSTM.pkl' 
if __name__ == '__main__': 
 args = ms_args_parser() 
 flag = 'ms' 
 train(args, LSTM_PATH, flag) 
 test(args, LSTM_PATH, flag) 
## args.py 
import argparse 
import torch 
def ms_args_parser(): 
 parser = argparse.ArgumentParser() 
 parser.add_argument('--epochs', type=int, default=30, help='input dimension') 
 parser.add_argument('--input_size', type=int, default=6, help='input dimension') 
 parser.add_argument('--output_size', type=int, default=21, help='output dimension') 
 parser.add_argument('--hidden_size', type=int, default=64, help='hidden size') 
 parser.add_argument('--num_layers', type=int, default=1, help='num layers') 
 parser.add_argument('--lr', type=float, default=0.05, help='learning rate') 
 parser.add_argument('--batch_size', type=int, default=1, help='batch size') 
 parser.add_argument('--optimizer', type=str, default='adam', help='type of optimizer') 
 parser.add_argument('--device', default=torch.device("cuda" if torch.cuda.is_available() else  
"cpu")) 
 parser.add_argument('--weight_decay', type=float, default=1e-4, help='weight decay') 
 parser.add_argument('--bidirectional', type=bool, default=False, help='LSTM direction') 
 args = parser.parse_args() 
 return args 
35## data_process.py 
import os 
import random 
import numpy as np 
import pandas as pd 
import torch 
from torch.utils.data import Dataset, DataLoader 
device = torch.device("cuda" if torch.cuda.is_available() else "cpu") 
def setup_seed(seed): 
 os.environ['PYTHONHASHSEED'] = str(seed) 
 torch.manual_seed(seed) 
 torch.cuda.manual_seed_all(seed) 
 np.random.seed(seed) 
 random.seed(seed) 
 torch.backends.cudnn.deterministic = True 
def load_data(): 
 path = os.path.dirname(os.path.realpath(__file__)) + '/data/q2(200).csv' 
 df = pd.read_csv(path, encoding='gbk') 
 columns = df.columns 
 df.fillna(df.mean(), inplace=True) 
 MAX = np.max(df[columns[0]]) 
 MIN = np.min(df[columns[0]]) 
 df[columns[0]] = (df[columns[0]] - MIN) / (MAX - MIN) 
 return df, MAX, MIN 
class MyDataset(Dataset): 
 def __init__(self, data): 
 self.data = data 
 def __getitem__(self, item): 
 return self.data[item] 
 def __len__(self): 
 return len(self.data) 
def nn_seq_ms(B): 
 print('data processing...') 
36 data, m, n = load_data() 
 load = data[data.columns[0]] 
 load = load.tolist() 
 data = data.values.tolist() 
 seq = [] 
 for i in range(len(data) - 50 - 21): 
 train_seq = [] 
 train_label = [] 
 for j in range(i, i + 50): 
 x = [load[j]] 
 for c in range(1, 6): 
 x.append(data[j][c]) 
 train_seq.append(x) 
 train_label.append(load[i + 50: i + 71]) 
 train_seq = torch.FloatTensor(train_seq) 
 train_label = torch.FloatTensor(train_label).view(-1) 
 seq.append((train_seq, train_label)) 
 Dtr = seq[0:int(len(seq))] 
 seq = [] 
 i = len(data) - 50 
 train_seq = [] 
 for j in range(i, i + 50): 
 x = [load[j]] 
 for c in range(1, 6): 
 x.append(data[j][c]) 
 train_seq.append(x) 
 train_seq = torch.FloatTensor(train_seq) 
 seq.append(train_seq) 
 Dte = seq[0:int(len(seq))] 
 # Dtr = seq[0:int(len(seq) * 0.7)] 
 # Dte = seq[int(len(seq) * 0.7):len(seq)] 
 # train_len = int(len(Dtr) / B) * B 
 # test_len = int(len(Dte) / B) * B 
 # Dtr, Dte = Dtr[:train_len], Dte[:test_len] 
 train = MyDataset(Dtr) 
 test = MyDataset(Dte) 
 Dtr = DataLoader(dataset=train, batch_size=B, shuffle=False, num_workers=0) 
 Dte = DataLoader(dataset=test, batch_size=B, shuffle=False, num_workers=0) 
37 return Dtr, Dte, m, n 
def get_mape(x, y): 
 """ 
 :param x:true 
 :param y:pred 
 :return:MAPE 
 """ 
 return np.mean(np.abs((x - y) / x)) 
## models.py 
from torch import nn 
import torch 
device = torch.device("cuda" if torch.cuda.is_available() else "cpu") 
class LSTM(nn.Module): 
 def __init__(self, input_size, hidden_size, num_layers, output_size, batch_size): 
 super().__init__() 
 self.input_size = input_size 
 self.hidden_size = hidden_size 
 self.num_layers = num_layers 
 self.output_size = output_size 
 self.num_directions = 1 
 self.batch_size = batch_size 
 self.lstm = nn.LSTM(self.input_size, self.hidden_size, self.num_layers,  
batch_first=True) 
 self.linear = nn.Linear(self.hidden_size, self.output_size) 
 def forward(self, input_seq): 
 h_0 = torch.randn(self.num_directions * self.num_layers, self.batch_size,  
self.hidden_size).to(device) 
 c_0 = torch.randn(self.num_directions * self.num_layers, self.batch_size,  
self.hidden_size).to(device) 
 # print(input_seq.size()) 
 seq_len = input_seq.shape[1] 
 # seq_len = 21 
 # input(batch_size, seq_len, input_size) 
 input_seq = input_seq.view(self.batch_size, seq_len, self.input_size) 
 # output(batch_size, seq_len, num_directions * hidden_size) 
 output, _ = self.lstm(input_seq, (h_0, c_0)) 
 # print('output.size=', output.size()) 
 # print(self.batch_size * seq_len, self.hidden_size) 
38 output = output.contiguous().view(self.batch_size * seq_len, self.hidden_size) # (5 *  
30, 64) 
 pred = self.linear(output) # pred() 
 # print('pred=', pred.shape) 
 pred = pred.view(self.batch_size, seq_len, -1) 
 pred = pred[:, -1, :] 
 return pred 
## util.py 
from itertools import chain 
import torch 
from scipy.interpolate import make_interp_spline 
from torch import nn 
import numpy as np 
import matplotlib.pyplot as plt 
from models import LSTM, BiLSTM 
from data_process import nn_seq, nn_seq_ms, nn_seq_mm, device, get_mape, setup_seed 
setup_seed(20) 
def train(args, path, flag): 
 if flag == 'us': 
 Dtr, Dte, m, n = nn_seq(B=args.batch_size) 
 elif flag == 'ms': 
 Dtr, Dte, m, n = nn_seq_ms(B=args.batch_size) 
 else: 
 Dtr, Dte, m, n = nn_seq_mm(B=args.batch_size, num=args.output_size) 
 input_size, hidden_size, num_layers = args.input_size, args.hidden_size, args.num_layers 
 output_size = args.output_size 
 if args.bidirectional: 
 model = BiLSTM(input_size, hidden_size, num_layers, output_size,  
batch_size=args.batch_size).to(device) 
 else: 
 model = LSTM(input_size, hidden_size, num_layers, output_size,  
batch_size=args.batch_size).to(device) 
 loss_function = nn.MSELoss().to(device) 
 if args.optimizer == 'adam': 
 optimizer = torch.optim.Adam(model.parameters(), lr=args.lr, 
 weight_decay=args.weight_decay) 
3940 
 else: 
 optimizer = torch.optim.SGD(model.parameters(), lr=args.lr, 
 momentum=0.9, weight_decay=args.weight_decay) 
 # training 
 loss = 0 
 for i in range(args.epochs): 
 cnt = 0 
 for (seq, label) in Dtr: 
 cnt += 1 
 seq = seq.to(device) 
 label = label.to(device) 
 y_pred = model(seq) 
 loss = loss_function(y_pred, label) 
 optimizer.zero_grad() 
 loss.backward() 
 optimizer.step() 
 # if cnt % 100 == 0: 
 # print('epoch', i, ':', cnt - 100, '~', cnt, loss.item()) 
 print('epoch', i, ':', loss.item()) 
 state = {'model': model.state_dict(), 'optimizer': optimizer.state_dict()} 
 torch.save(state, path) 
def test(args, path, flag): 
 if flag == 'us': 
 Dtr, Dte, m, n = nn_seq(B=args.batch_size) 
 elif flag == 'ms': 
 Dtr, Dte, m, n = nn_seq_ms(B=args.batch_size) 
 else: 
 Dtr, Dte, m, n = nn_seq_mm(B=args.batch_size, num=args.output_size) 
 pred = [] 
 y = [] 
 print('loading model...') 
 input_size, hidden_size, num_layers = args.input_size, args.hidden_size, args.num_layers 
 output_size = args.output_size 
 if args.bidirectional: 
 model = BiLSTM(input_size, hidden_size, num_layers, output_size,  
batch_size=args.batch_size).to(device) 
 else: 
 model = 
LSTM(input_size, hidden_size, num_layers, output_size,  
batch_size=args.batch_size).to(device) 
 # model = LSTM(input_size, hidden_size, num_layers, output_size,  
batch_size=args.batch_size).to(device) model.load_state_dict(torch.load(path)['model']) 
 model.eval() 
 print('predicting...') 
 for (seq, target) in Dte: 
 target = list(chain.from_iterable(target.data.tolist())) 
 y.extend(target) 
 seq = seq.to(device) 
 with torch.no_grad(): 
 y_pred = model(seq) 
 y_pred = list(chain.from_iterable(y_pred.data.tolist())) 
 pred.extend(y_pred) 
 for seq in Dte: 
 seq = seq.to(device) 
 with torch.no_grad(): 
 y_pred = model(seq) 
 y_pred = list(chain.from_iterable(y_pred.data.tolist())) 
 pred.extend(y_pred) 
 y, pred = np.array([y]), np.array([pred]) 
 # pred =np.array([pred]) 
 y = (m - n) * y + n 
 pred = (m - n) * pred + n 
 print(pred) 
 print('mape:', get_mape(y, pred)) 
 # plot 
 x = [i for i in range(len(y[0]))] 
 # x_smooth = np.linspace(np.min(x), np.max(x), 900) 
 # y_smooth = make_interp_spline(x, y)(x_smooth) 
 plt.plot(x, y[0][0:len(y[0])], c='green', marker='*', ms=1, alpha=0.75, label='true') 
 # y_smooth = make_interp_spline(x, pred.T[150:300])(x_smooth) 
 plt.plot(x, pred[0][0:len(pred[0])], c='red', marker='o', ms=1, alpha=0.75, label='pred') 
 plt.grid(axis='y') 
 plt.legend() 
 plt.show() 
第三问 
aSOC<-read.table(file="C:/Users/DELL/Desktop/NG-G17.csv",sep=",",header=T) 
xSOC<-ts(aSOC$SOC,start=2012) 
auto.arima(xSOC) 
SOC.fit<-arima(xSOC,order=c(0,0,0)) 
SOC.fore<-forecast(SOC.fit,h=1) 
SOC.fore 
41aSIC<-read.table(file="C:/Users/DELL/Desktop/NG-G17.csv",sep=",",header=T) 
xSIC<-ts(aSIC$SIC,start=2012) 
auto.arima(xSIC) 
xSIC.fit<-arima(xSIC,order=c(0,0,0)) 
xSIC.fore<-forecast(xSIC.fit,h=1) 
xSIC.fore 
aSTC<-read.table(file="C:/Users/DELL/Desktop/NG-G17.csv",sep=",",header=T) 
xSTC<-ts(aSTC$STC,start=2012) 
auto.arima(xSTC) 
xSTC.fit<-arima(xSTC,order=c(0,0,0)) 
xSTC.fore<-forecast(xSTC.fit,h=1) 
xSTC.fore 
aN<-read.table(file="C:/Users/DELL/Desktop/NG-G17.csv",sep=",",header=T) 
xN<-ts(aN$N,start=2012) 
auto.arima(xN) 
xN.fit<-arima(xN,order=c(0,0,0)) 
xN.fore<-forecast(xN.fit,h=1) 
xN.fore 
aCN<-read.table(file="C:/Users/DELL/Desktop/NG-G17.csv",sep=",",header=T) 
xCN<-ts(aCN$CN,start=2012) 
auto.arima(xCN) 
xCN.fit<-arima(xCN,order=c(0,0,0)) 
xCN.fore<-forecast(xCN.fit,h=1) 
xCN.fore 
第四问 
a<-read.table(file="C:/Users/DELL/Desktop/主成分分析.csv",sep=",",header=T) 
x1<-ts(a$x1) 
x2<-ts(a$x2) 
x3<-ts(a$x3) 
x4<-ts(a$x4) 
x5<-ts(a$x5) 
x6<-ts(a$x6) 
X<-data.frame(x1,x2,x3,x4,x5,x6) 
test.pr<-princomp(X,cor=FALSE) 
summary(test.pr,loadings=TRUE) #做主成分分析 
plot(test.pr,type="line") #画主成分碎石图 
a<-read.table(file="C:/Users/DELL/Desktop/NG.csv",sep=",",header=T) 
x1<-ts(a$sd10) 
x2<-ts(a$sd40) 
42x3<-ts(a$sd100) 
x4<-ts(a$sd200) 
x5<-ts(a$CN) 
X<-data.frame(x1,x2,x3,x4,x5) 
test.pr<-princomp(X) 
summary(test.pr,loadings=TRUE) #做主成分分析 
plot(test.pr,type="line") #画主成分碎石图 
第六问 
aSOC<-read.table(file="C:/Users/DELL/Desktop/6.csv",sep=",",header=T) 
xSOC<-ts(aSOC$SOC) 
auto.arima(xSOC) 
SOC.fit<-arima(xSOC,order=c(0,1,0)) 
SOC.fore<-forecast(SOC.fit,h=2) 
SOC.fore 
aSIC<-read.table(file="C:/Users/DELL/Desktop/6.csv",sep=",",header=T) 
xSIC<-ts(aSIC$SIC) 
auto.arima(xSIC) 
xSIC.fit<-arima(xSIC,order=c(0,1,0)) 
xSIC.fore<-forecast(xSIC.fit,h=2) 
xSIC.fore 
aSTC<-read.table(file="C:/Users/DELL/Desktop/6.csv",sep=",",header=T) 
xSTC<-ts(aSTC$STC) 
auto.arima(xSTC) 
xSTC.fit<-arima(xSTC,order=c(0,0,0)) 
xSTC.fore<-forecast(xSTC.fit,h=2) 
xSTC.fore 
aN<-read.table(file="C:/Users/DELL/Desktop/6.csv",sep=",",header=T) 
xN<-ts(aN$N) 
auto.arima(xN) 
xN.fit<-arima(xN,order=c(0,1,0)) 
xN.fore<-forecast(xN.fit,h=2) 
xN.fore 
aCN<-read.table(file="C:/Users/DELL/Desktop/6.csv",sep=",",header=T) 
xCN<-ts(aCN$CN) 
auto.arima(xCN) 
xCN.fit<-arima(xCN,order=c(0,1,0)) 
xCN.fore<-forecast(xCN.fit,h=2) 
xCN.fore 
43aNDVI<-read.table(file="C:/Users/DELL/Desktop/6.csv",sep=",",header=T) 
xNDVI<-ts(aNDVI$NDVI) 
auto.arima(xNDVI) 
xNDVI.fit<-arima(xNDVI,order=c(2,0,2)) 
xNDVI.fore<-forecast(xNDVI.fit,h=17) 
xNDVI.fore