数学建模，8月学习感想

长文预警！！！    

Matlab入门

   Matlab可谓是数学建模中的热门语言，其开发模式可以分为命令行模式、脚本模式和面向对象模式,其基本语法同大多数编程语言相似。  可参照   Matlab语法教程

   由于绝大部分的数学建模问题都可以通过Matlab编程后求解，因此，了解Matlab的使用至关重要。
 

常用的操作指令

操作指令	功能
clc	清除指令窗口中显示的内容
clear	清除MATLAB工作空间中保存的变量
close all	关闭所有打开的图形窗口
clf	清楚图形窗的内容
edit	打开m文件编辑器
disp	显示变量的内容

数据类型

Matlab常用的数据类型有数值型（主要应用整型和浮点型）、字符型、函数句柄、元胞数组和结构体数组，此外也有映射表（Map）、表（Table）等高级数据类型。
 
关于常用数据类型部分，这里有一篇讲的详细的文章。
Matlab数据类型

plus：学会利用Matlab封装好的各类工具箱,也能简化我们的问题求解

  

建模流程

   分析问题去建立可求解的模型往往是建模中最为难产的过程 （脑袋发懵属于典型症状） ，虽然论文的撰写也需要在广泛阅读名家论文基础上的苦心经营。

数学建模的基本流程如下：

  

分析问题

如何分析问题体现了我们建模的思路，这跟做数学应用题的感觉差不多。由于是个脑力活，所以需要经历头脑风暴，在广泛查阅各项资料的基础上将问题考虑得 又全面又高效 最好。分析问题往往是建立模型的基础，但建立的模型给问题的分析负反馈时，也总能得到新的思路，所以注意相互的印证。
  
如“将数据保存为表格，再将表格抽象为矩阵”这样的一种从数据分析入手的正向思路； 也有像借鉴动态规划的分阶段思想，先将大问题分解为小问题，通过将全过程划分为一个个相互联系的局部规划，最终达到整体最优的分解思路； 还有如在优化模型分析中抽象问题为决策变量（寻求的决策是什么）、目标函数（需要优化的目标是什么）和约束条件（决策面临哪些限制条件），抽象所给条件后，寻找与已有经典模型间距离的类比思路； 也有从最后我们所要给出理想化模型出发，假设我们是要建立一个离散模型，则我们需要做的是在数据离散化后再通过路径规划来分析的一种反向思路。
  

建立模型&求解模型

常见的数学模型要建立在 可求解和可说明 的基础上。可求解，是指它能通过你的求解工具算出答案，可说明，是指它能体现你的分析思路，并能在论文部分说明。
  
常见的数据建模方法可以分为数据建模技术、优化技术、连续模型求解、评价模型求解和机理建模方法¹。其中数据建模技术和优化技术种类繁多，机理建模方法往往需要另辟蹊径，都值得倍加留意。
  

数据建模技术

  
在大数据时代，数据分析已然成为热门。一张表，一个excel，其实已经能满足大部分人的日常数据分析需求。但面临数学建模时，你可能需要更高级的代码支持，此时，Matlab和Python就都是你的可选之才。  虽然工具还有很多，如SPSS
  
对于Python，主要通过numpy、pandas、scikit-learn等科学计算库实现对数据的分析，又主要通过matplotlib等绘图库实现数据的可视化。
  
常见的数据分析包括数据准备（对数据进行标准化、归一化或二值化处理），降维（PCA、随机投影、NMF等，目的是将原高维空间中的数据点映射到低纬度空间，降低分析量），聚类（K-Means，层次聚类，GMM等），分类（Logistic回归，支持向量机，朴素贝叶斯等），回归（最小二乘回归，决策树回归，岭回归等）。值得注意的是，数据分析往往是机器学习的基础，因此，掌握基本的数据建模技术，将不止步于数据建模领域。 比如开启自主学习的大门
  

下面以一个最小二乘法的的实例来说明，实现的是数据的线性回归²。
首先是Matlab的代码实现部分：

%% 输入数据，可替换
clc, clear all, close all
x=[23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40];
y=[41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0];	
%% 采用最小二乘回归
% 绘制散点图，判断是否具有线性关系
figure
plot(x,y,'r*')                         %作散点图
xlabel('x（职工工资总额）','fontsize', 12)           %横坐标名
ylabel('y（商品零售总额）', 'fontsize',12)           %纵坐标名
set(gca,'linewidth',2);
% 采用最小二乘拟合
Lxx=sum((x-mean(x)).^2);
Lxy=sum((x-mean(x)).*(y-mean(y)));
b1=Lxy/Lxx;
b0=mean(y)-b1*mean(x);
y1=b1*x+b0;
hold on
plot(x, y1,'linewidth',2);
fprintf('b1= %d   bo=%d \n',b1,b0);
fprintf('y1=%d*x+(%d) \n',b1,b0);

拟合结果：

b1= 2.799121e+00   bo=-2.354935e+01 
y1=2.799121e+00*x+(-2.354935e+01) 

然后是比较类似的Python代码部分，对比可以发现：Matlab最大的特点是自身集成好了我们运算所用到的库，而在Python中则需要导入相应的库后，才能使用某方法。

#-*- coding:utf-8 -*-

##最小二乘法
import numpy as np   ##科学计算库
import matplotlib.pyplot as plt  ##绘图库
from scipy.optimize import leastsq  ##引入最小二乘法算法
import pylab as mpl

'''
     设置中文显示字体，因为matplotlib默认不支持中文
'''
mpl.rcParams['font.sans-serif'] = ['FangSong']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题


'''
     设置样本数据，真实数据需要在这里处理
'''
##样本数据(Xi,Yi)，需要转换成数组(列表)形式
Xi=np.array([23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40])
Yi=np.array([41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0])

'''
    设定拟合函数和偏差函数
    函数的形状确定过程：
    1.先画样本图像
    2.根据样本图像大致形状确定函数形式(直线、抛物线、正弦余弦等)
'''

##需要拟合的函数func :指定函数的形状
def func(p,x):
    k,b=p
    return k*x+b

##偏差函数：x,y都是列表:这里的x,y更上面的Xi,Yi中是一一对应的
def error(p,x,y):
    return func(p,x)-y

#k,b的初始值，可以任意设定,经过几次试验，发现p0的值会影响cost的值：Para[1]
p0=[1,20]

#把error函数中除了p0以外的参数打包到args中(使用要求)
Para=leastsq(error,p0,args=(Xi,Yi))

#读取结果
k,b=Para[0]
print("k=",k,"b=",b)
print("cost："+str(Para[1]))
print("求解的拟合直线为:")
print("y="+str(round(k,2))+"x+("+str(round(b,2))+")")

'''
   绘图，看拟合效果.
'''

#画样本点
plt.figure(figsize=(8,6)) ##指定图像比例： 8：6
plt.scatter(Xi,Yi,color="red",label="样本数据",linewidth=2)

#画拟合直线
x=np.linspace(20,80,300) ##在20-80直接画300个连续点
y=k*x+b ##函数式
plt.plot(x,y,color="red",label="拟合直线",linewidth=2)
plt.legend(loc='lower right') #绘制图例
plt.xlabel('x（职工工资总额）')#横坐标
plt.ylabel('y（商品零售总额）')#纵坐标
plt.show()

拟合结果：

k= 2.7991214484021274 b= -23.54934679295249
cost：1
求解的拟合直线为:
y=2.8x+(-23.55)

（结果相差无几，只是绘图风格略有差异）
 
 

优化技术

当说到优化技术，其实我们应该先提到运筹学。在百度百科上，显示运筹学的主要内容有：

1. 规划论：包括线性规划、非线性规划、整数规划和动态规划 （即数学规划课程的一般内容） 。线性规划  以单纯形法的提出而闻名，本质上是求解线性方程组的问题。由于空间与解集间存在对应性，使得通过行列式和矩阵的迭代就能明确规划是否可解，虽然当约束条件数过多后，将不得不借助计算机来加速计算进程。线性规划与非线性规划的主要区别在于约束条件和目标函数是否呈现为线性关系。对于 非线性规划 的求解， 其会更多地依赖于传统微分学的知识（二次规划）和搜索算法（牛顿法），而如蚁群算法、模拟退火算法、遗传算法等启发式算法其实与这部分的内容兼容。整数规划  的最大特点来源于整数，这使得可行域由连续空间变为离散的点，分支定界法和割平面法是其中两种比较经典的算法，而0-1问题则属于整数规划中的特殊问题。动态规划  指的是在多阶段过程中的决策序列问题，背包问题和旅行商问题就属于其中的两个典型例子，“阶段”、“状态”、“无后效性”、“决策”、“策略”是动态规划中重要的术语概念。值得注意的是，许多优化算法往往并不能给出全局的最优解，而只能取到一个尽可能好的解，这同贪心算法的逻辑一致，毕竟资源有限的情况下我们并没那么多时间穷举所有的可能性，甚至对于计算机而言。
    
2. 库存论：库存论是一种研究物质最优存储及存储控制的理论，物质存储是工业生产和经济运转的必然现象。如果物质存储过多，则会占用大量仓储空间，增加保管费用，使物质过时报废从而造成经济损失；如果存储过少，则会因失去销售时机而减少利润，或因原料短缺而造成停产。虽然物品的数量通常以离散的形式出现，但在求解时我们往往考虑转化成连续模型便于分析。
    
3. 图论：在解决工程设计和管理决策的最优化问题中发挥作用（可以把动态规划中的状态想象成一个节点），经典问题如哥德斯堡七桥问题，经典算法如Dijkstra算法和Floyd算法。
    
4. 排队论：又称随机服务系统理论，主要研究各种系统的排队长度，排队的等待时间和所提供的服务等各种参数，给出系统随机聚散现象的解释。在生活中往往具有广泛的应用场景，如水库水量的调节、生产流水线的安排、铁路分成场的调度、电网的设计等。
    
5. 对策论：又称博弈论，是主要研究双方冲突、制胜对策的理论。经典问题如田忌赛马。经典理论如追求“你好我好大家好”的纳什均衡。

一个博弈过程中，无论对方的策略选择如何，当事人一方都会选择某个确定的策略，则该策略被称作支配性策略。如果任意一位参与者在其他所有参与者的策略确定的情况下，其选择的策略是最优的，那么这个组合就被定义为纳什平衡。

对于数学规划问题，大部分都能够用Lingo方便解决。而对于常见的优化算法，Matlab和Python其实都能提出自己的解决方案。
 
 
如用Lingo求解一个线性规划：

max = 4*x1 + 2*x2 + 5*x3;
x1 + 2*x2 + x3 <=20;
2*x1 + 3*x2 + x3 >= 14;
x1 + 2*x2 + x3 = 10;
@free(x3);

 
结果：
 
 
再如用matlab实现蚁群算法：

%% 蚁群算法及Matlab实现——TSP问题

%% 数据准备
% 清空环境变量
clear all
clc

% 程序运行计时开始
t0 = clock;

%导入数据
citys=xlsread('Chap9_citys_data.xlsx', 'B2:C53');

%% 计算城市间相互距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));
        else
            D(i,j) = 1e-4;      %设定的对角矩阵修正值
        end
    end    
end

%% 初始化参数
m = 75;                              % 蚂蚁数量
alpha = 1;                           % 信息素重要程度因子
beta = 5;                            % 启发函数重要程度因子
vol = 0.2;                           % 信息素挥发(volatilization)因子
Q = 10;                              % 常系数
Heu_F = 1./D;                        % 启发函数(heuristic function)
Tau = ones(n,n);                     % 信息素矩阵
Table = zeros(m,n);                  % 路径记录表
iter = 1;                            % 迭代次数初值
iter_max = 100;                      % 最大迭代次数 
Route_best = zeros(iter_max,n);      % 各代最佳路径       
Length_best = zeros(iter_max,1);     % 各代最佳路径的长度  
Length_ave = zeros(iter_max,1);      % 各代路径的平均长度  
Limit_iter = 0;                      % 程序收敛时迭代次数

%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
      start = zeros(m,1);
      for i = 1:m
          temp = randperm(n);
          start(i) = temp(1);
      end
      Table(:,1) = start; 
      % 构建解空间
      citys_index = 1:n;
      % 逐个蚂蚁路径选择
      for i = 1:m
          % 逐个城市路径选择
         for j = 2:n
             tabu = Table(i,1:(j - 1));           % 已访问的城市集合(禁忌表)
             allow_index = ~ismember(citys_index,tabu);    % 参加说明1（程序底部）
             allow = citys_index(allow_index);  % 待访问的城市集合
             P = allow;
             % 计算城市间转移概率
             for k = 1:length(allow)
                 P(k) = Tau(tabu(end),allow(k))^alpha * Heu_F(tabu(end),allow(k))^beta;
             end
             P = P/sum(P);
             % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);     %参加说明2(程序底部)
            target_index = find(Pc >= rand); 
            target = allow(target_index(1));
            Table(i,j) = target;
         end
      end
      % 计算各个蚂蚁的路径距离
      Length = zeros(m,1);
      for i = 1:m
          Route = Table(i,:);
          for j = 1:(n - 1)
              Length(i) = Length(i) + D(Route(j),Route(j + 1));
          end
          Length(i) = Length(i) + D(Route(n),Route(1));
      end
      % 计算最短路径距离及平均距离
      if iter == 1
          [min_Length,min_index] = min(Length);
          Length_best(iter) = min_Length;  
          Length_ave(iter) = mean(Length);
          Route_best(iter,:) = Table(min_index,:);
          Limit_iter = 1; 
          
      else
          [min_Length,min_index] = min(Length);
          Length_best(iter) = min(Length_best(iter - 1),min_Length);
          Length_ave(iter) = mean(Length);
          if Length_best(iter) == min_Length
              Route_best(iter,:) = Table(min_index,:);
              Limit_iter = iter; 
          else
              Route_best(iter,:) = Route_best((iter-1),:);
          end
      end
      % 更新信息素
      Delta_Tau = zeros(n,n);
      % 逐个蚂蚁计算
      for i = 1:m
          % 逐个城市计算
          for j = 1:(n - 1)
              Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);
          end
          Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);
      end
      Tau = (1-vol) * Tau + Delta_Tau;
    % 迭代次数加1，清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end

%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost=etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost) '秒']);

%% 绘图
figure(1)
plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],...  %三点省略符为Matlab续行符
     [citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-');
grid on
for i = 1:size(citys,1)
    text(citys(i,1),citys(i,2),['   ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),'       起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),'       终点');
xlabel('城市位置横坐标')
ylabel('城市位置纵坐标')
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')'])
figure(2)
plot(1:iter_max,Length_best,'b')
legend('最短距离')
xlabel('迭代次数')
ylabel('距离')
title('算法收敛轨迹')
%--------------------------------------------------------------------------
%% 程序解释或说明
% 1. ismember函数判断一个变量中的元素是否在另一个变量中出现，返回0-1矩阵；
% 2. cumsum函数用于求变量中累加元素的和，如A=[1, 2, 3, 4, 5], 那么cumsum(A)=[1, 3, 6, 10, 15]。

结果如下：

最短距离:7795.6696
最短路径:21  23  20  50  16  29  46  44  34  35  36  39  40  38  37  48  24   5  15   6   4  25  12  28  27  26  47  13  14  52  11  51  33  43  10   9   8  41  19  45  32  49   1  22  31  18   3  17  42   7   2  30  21
收敛迭代次数:30
程序执行时间:5.251秒

 
 
再如用python实现模拟退火算法，可参照这篇文章：模拟退火算法解释
解决如下极值问题：$$maxF=5x^2-3.4x^4+\frac{x^6}{5}+xy-4y^3+2y^4$$

import math
from random import random
import matplotlib.pyplot as plt


def func(x, y):  # 函数优化问题
    res = 5 * x ** 2 - 3.4 * x ** 4 + x ** 6 / 5 + x * y - 4 * y ** 3 + 2 * y ** 4
    return res

#定义了SA类
# x为公式里的x1,y为公式里面的x2
class SA:
    def __init__(self, func, iter=200, T0=100, Tf=0.01, alpha=0.99):
        self.func = func
        self.iter = iter  # 内循环迭代次数,即为L =200
        self.alpha = alpha  # 降温系数，alpha=0.99
        self.T0 = T0  # 初始温度T0为100
        self.Tf = Tf  # 温度终值Tf为0.01
        self.T = T0  # 当前温度
        self.x = [random() * 11 - 5 for i in range(iter)]  # 随机生成200个x的值
        self.y = [random() * 11 - 5 for i in range(iter)]  # 随机生成200个y的值
        self.most_best = []
        self.history = {'f': [], 'T': []}

    def generate_new(self, x, y):  # 扰动产生新解的过程
        while True:
            x_new = x + self.T * (random() - random())
            y_new = y + self.T * (random() - random())
            if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
                break  # 重复得到新解，直到产生的新解满足约束条件
        return x_new, y_new

    def Metrospolis(self, f, f_new):  # Metropolis准则
        if f_new <= f:
            return 1
        else:
            p = math.exp((f - f_new) / self.T)
            if random() < p:
                return 1
            else:
                return 0

    def best(self):  # 获取最优目标函数值
        f_list = []  # f_list数组保存每次迭代之后的值
        for i in range(self.iter):
            f = self.func(self.x[i], self.y[i])
            f_list.append(f)
        f_best = min(f_list)

        idx = f_list.index(f_best)
        return f_best, idx  # f_best,idx分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下标

    def run(self):
        count = 0
        # 外循环迭代，当前温度小于终止温度的阈值
        while self.T > self.Tf:

            # 内循环迭代100次
            for i in range(self.iter):
                f = self.func(self.x[i], self.y[i])  # f为迭代一次后的值
                x_new, y_new = self.generate_new(self.x[i], self.y[i])  # 产生新解
                f_new = self.func(x_new, y_new)  # 产生新值
                if self.Metrospolis(f, f_new):  # 判断是否接受新值
                    self.x[i] = x_new  # 如果接受新值，则把新值的x,y存入x数组和y数组
                    self.y[i] = y_new
            # 迭代L次记录在该温度下最优解
            ft, _ = self.best()
            self.history['f'].append(ft)
            self.history['T'].append(self.T)
            # 温度按照一定的比例下降（冷却）
            self.T = self.T * self.alpha
            count += 1

            # 得到最优解
        f_best, idx = self.best()
        print(f"  解F={f_best}, x={self.x[idx]}, y={self.y[idx]}")


sa = SA(func)
sa.run()

plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()

结果如下：

解  F=-99.33763868366601, x=-3.2523049021920905, y=1.6481055790732624

 
 

连续模型求解

 
所谓连续模型，主要是在求解微分方程（常微分方程和偏微分方程），举个例子:

已知matlab的slovepde函数可以求解如下形式的方程：
$$m\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot (c\nabla u）+au=f$$

于是只要设置好参数后，就可以求解这个二阶波的方程：
$$2\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot \nabla u=0$$

%% 设置参数
c = 1;
a = 0;
f = 0;
m = 2;

%% 定义波的空间位置
numberOfPDE = 1;
model = createpde(numberOfPDE);
geometryFromEdges(model,@squareg);
pdegplot(model,'EdgeLabels','on'); 
ylim([-1.1 1.1]);
axis equal
title 'Geometry With Edge Labels Displayed';
xlabel x
ylabel y

%% 定义微分方程模型的系数和边界条件
specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);
applyBoundaryCondition(model,'dirichlet','Edge',[2,4],'u',0);
applyBoundaryCondition(model,'neumann','Edge',([1 3]),'g',0);

%% 定义该问题的有限元网格
generateMesh(model);
figure
pdemesh(model);
ylim([-1.1 1.1]);
axis equal
xlabel x
ylabel y

%% 定义初始条件
u0 = @(location) atan(cos(pi/2*location.x));
ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));
setInitialConditions(model,u0,ut0);

%% 方程求解
n = 29;
tlist = linspace(0,5,n);
model.SolverOptions.ReportStatistics ='on';
result = solvepde(model,tlist);
u = result.NodalSolution;

%% 模型的数值仿真
figure
umax = max(max(u));
umin = min(min(u));
for i = 1:n
    pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ZStyle','continuous',...
                  'Mesh','off','XYGrid','on','ColorBar','off');
    axis([-1 1 -1 1 umin umax]); 
    caxis([umin umax]);
    xlabel x
    ylabel y
    zlabel u
    M(i) = getframe;
end

运行后，就可以得到波的数值仿真效果图：

 
 

评价模型求解

 
评价模型多用于在多个选项间做出评价分析，关键在于权重矩阵的计算，常用的两种方法是线性加权法和层次分析法。
 
线性加权法的关键是构建评价的指标体系，因为权重较大指标值对综合指标作用较大，所以需综合考虑评价指标的权重值。
 
层次分析法（AHP）是一种基础的算法，在评价评判、资源分配和决策和优化问题中都有应用。matlab主要实现将评判矩阵转化为因素的权重矩阵的过程，评判矩阵则依据个人在两两比较中的权重分配。
$$a_{ij}表示与j指标相比，i指标的重要程度，值越大越重要，且满足a_{ij}\cdot a_{ji}=1$$
$$评判矩阵\left(\begin{array}{ccc}1& 2& 6\\ & & \\ \frac{1}{2}& 1& 4\\ & & \\ \frac{1}{6}& \frac{1}{4}& 1\end{array}\right)$$
给出层次分析法的matlab代码示例。

% 数据读入
clc
clear all
A=[1 2 6; 1/2 1 4; 1/6 1/4 1];% 评判矩阵
% 一致性检验和权向量计算
[n,n]=size(A);
[v,d]=eig(A);
r=d(1,1);
CI=(r-n)/(n-1);
RI=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.52 1.54 1.56 1.58 1.59];
CR=CI/RI(n);
if  CR<0.10
    CR_Result='通过';
   else
    CR_Result='不通过';   
end

% 权向量计算
w=v(:,1)/sum(v(:,1));
w=w';

%% 结果输出
disp('该判断矩阵权向量计算报告：');
disp(['一致性指标:' num2str(CI)]);
disp(['一致性比例:' num2str(CR)]);
disp(['一致性检验结果:' CR_Result]);
disp(['特征值:' num2str(r)]);
disp(['权向量:' num2str(w)]);

结果如下，可见1的权重最大，2、3其次：

一致性指标:0.0046014
一致性比例:0.0079334
一致性检验结果:通过
特征值:3.0092
权向量:0.58763     0.32339    0.088983

 

机理建模方法

 
机理建模往往具有明确的物理意义，模型参数易于调整，在遭遇一些非典型的数学建模问题（非数据、优化、连续、评价）时，适应性良好。缺点在于数理化的过程可能较为痛苦。虽然具有一定的思维难度，但掌握基本的机理建模思路和方法应该说是大势所趋。
 

机理建模，是根据系统的机理（如物理或化学的变化规律）建立系统模型的过程。首先，根据建模对象的应用场合和模型的使用目的进行合理的假设；随后，根据系统的内在机理建立数学方程，并比较过程变量数与独立方程数来进行自由度分析，以保证模型有解；最后，进行模型简化与验证。

 
机理建模常见的有两类，一类是推导法机理建模（基于物理过程），类似于微分方程建模，常用于动力学的建模过程，比如化学中反应动力学，还有各种场的方程，比如压力场、热场方程等；一类是包含一个或几个类别对象的复杂系统问题，常通过元胞自动机—仿真法来进行机理建模。
 

元胞自动机(cellular automata，CA) 是一种时间、空间、状态都离散，空间相互作用和时间因果关系为局部的网格动力学模型，具有模拟复杂系统时空演化过程的能力。

CA的基本思想是：CA会定义一个网格，网格上的每个点代表一个有限数量的状态中的细胞，过渡规则同时应用到每一个细胞。典型的转化规则是依赖细胞和它近邻的状态。或者可参照这篇文章入门：   元胞自动机入门
 
 
在实际建模中，要灵活定义元胞（事件的基本组件）、更新规则（元胞间的相互作用）和系统输出（以怎样的形式反映事件效果）。
 
如给出一个反映打车过程的元胞自动机matlab实现：

clc, clf,clear 

% 界面设计（环境的定义）
plotbutton=uicontrol('style','pushbutton',...
   'string','Run', ...
   'fontsize',12, ...
   'position',[100,400,50,20], ...
   'callback', 'run=1;');

% 定义 stop button
erasebutton=uicontrol('style','pushbutton',...
   'string','Stop', ...
   'fontsize',12, ...
   'position',[200,400,50,20], ...
   'callback','freeze=1;');

% 定义 Quit button
quitbutton=uicontrol('style','pushbutton',...
   'string','Quit', ...
   'fontsize',12, ...
   'position',[300,400,50,20], ...
   'callback','stop=1;close;');

number = uicontrol('style','text', ...
    'string','1', ...
   'fontsize',12, ...
   'position',[20,400,50,20]);

%  元胞自动机的设置
n=128;
z = zeros(n,n);
cells = z;
sum = z;
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;
cells = (rand(n,n))<.5 ;
imh = image(cat(3,cells,z,z));
axis equal
axis tight

% 元胞索引更新的定义
x = 2:n-1;
y = 2:n-1;

% 元胞更新的规则定义
stop= 0; % wait for a quit button push
run = 0; % wait for a draw 
freeze = 0; % wait for a freeze
while (stop==0) 
    if (run==1)
        %nearest neighbor sum
        sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
            cells(x-1, y) + cells(x+1,y) + ...
            cells(x-1,y-1) + cells(x-1,y+1) + ...
            cells(3:n,y-1) + cells(x+1,y+1);
        % The CA rule
        cells = (sum==3) | (sum==2 & cells);       
        % draw the new image
        set(imh, 'cdata', cat(3,cells,z,z) )
        % update the step number diaplay
        stepnumber = 1 + str2num(get(number,'string'));
        set(number,'string',num2str(stepnumber))
    end
    if (freeze==1)
        run = 0;
        freeze = 0;
    end
    drawnow  %need this in the loop for controls to work
end

初始状态，杂乱无章,：

在该规则下的运行状态，零散响应：

修改规则后的运行状态，响应规模变大：

 
  小结一下，数学建模的方法其实远不止于本文中所展示的那些，毕竟涉及到广泛的数学和计算机学知识后，知识已成为汪洋，按需学习才会是高效的方式。
 

撰写论文

  俗话说，巧妇难为无米之炊，当然是要建好了模型才能写论文。所以友情提示，写论文要先理清思路。也可参照人家讲的  清风建模写作
 
论文的一般结构如下：

   

前情提示

    这部分指的是标题、摘要和关键词部分，即论文给人的第一印象。往往题好一半文，因此一个能直指要害地说明所采取的建模方法和解决的数学问题的标题为人所欣赏。关键词既可以是建模过程中所采取的重要方法，也可以是描述该问题过程中的重要概念。
   摘要在前情提示中处于核心位置。一个好的摘要，要在能够体现模型的数学归类、建模思路、算法思路、建模特点和主要结果的基础上，言简意赅并完整独立，成为本篇文章的精华。因此往往需在文章的最后提炼，并在细细打磨后出品。
   

正文部分

   这里是对问题详细展开的部分，要力争做到内容翔实，不重不漏，同时尽量说得通俗晓畅，毕竟读者作者都是人，如果连自己都解释不清楚的东西，那就不要班门弄斧了。

1. 问题分析：描述问题背景、建模目的、搜集文献和数据的大致内容，注意要用自己的语言描述，不要照抄材料。  查重预警

2. 模型假设：根据题目中的条件和要求作出假设，同时关注模型自身的假设条件，如一些回归模型会对数据的分布情况提出要求，这就需要预先验证。关键性假设不能缺少，符号使用要简洁通用。当然，过于理想化的假设在降低难度的同时也会降低模型的适用性，“又简洁又高效”才是我们追求的标准。

3. 模型建立：首先给出基本模型，包括相应的数学公式和方案。如果模型有可以简化的地方，要说明简化的思路并给出简化模型。模型当然不是越复杂越好 （复杂≠高大上） ，而要实用有效，以解决问题为原则。能用简单通俗方法解决的，就不用复杂晦涩方法。鼓励一切基于切实有效的创新，提升科学性和推广性，但不可摆花架子。语言表达要客观明确，忌生造术语，保证文章的逻辑性和严谨性。

4. 模型求解：建立的数学命题要符合逻辑规范，采用的计算方法要点明基本原理，计算过程要抓住关键步骤，并尽量算出合理的数值结果。  实不相瞒，我们都是要看数据说话的

5. 模型检验：对数值结果或模拟结果进行必要的检验，并能解释相应的意义。对题目中的问题作出回答，并能说明方案的合理性。结果的表示要清晰集中，推荐使用图形图表等形式展示思路和结果。   一图胜千言

6. 模型评价：优点突出，缺点不回避。可通过灵敏度分析和误差分析说明模型的适应性。改变原题要求和推广模型应用，也可以在此处点明，但同样需要客观实际，不可天马行空。

7. 参考文献：参考文献一般基于科学范式，在正文中引用时也需点明。

被引用的文献为期刊论文的单篇文献时，著录格式为：“顺序号 作者 题号 刊名， 出版年， 卷号（期号），引文所在的起止页码”。
例：[1]方新贵，王敏.关于包装{(p,p-1)(p,p)}图对和Slater问题[J].系统科学与数学.1989;(9):133—137
 
被引用的文献为图书、科技报告等整本文献时，著录格式为：“顺序号 作者 文献书名 版本（第一版本不标注） 出版地址，出版者，出版年”。
例：[2]姜启源.数学模型[M].北京：高等教育出版社，2003.7—9

 

配套说明

   在附录部分，可以列出详细的数据表格和求解程序，保持逻辑清晰，排版美观。这可以算是图、表和代码的“重灾区”。

   
   

排版工具

   排版工具常用的为两种，Word和LaTex。其中Word “所见即所得”，大部分人在日常办公中早已上过手； LaTex“内容与格式分离”，通过宏包来编排文档格式，具有劝退门槛，虽然专业性应该更胜一筹。注意编译时采用XeLaTeX以支持中文。

   不过对于参加比赛的同学，往往有许多模版可供参考。因此掌握基本方法，学会构建与重组就好。此外，对于word排版，可以参照侯捷的《word排版的艺术》； 对于LaTex排版，可以参照《Ishort-cn》和刘海洋的《LaTex入门》。 比如参照这个LaTex的模板： LaTex国赛模板

   
当你输完这些代码：

\documentclass[withoutpreface,bwprint]{cumcmthesis}
\title{论文题目}
\begin{document}
\maketitle
\begin{abstract}
摘要的具体内容。
\keywords{关键词1\quad 关键词2\quad 关键词3}
\end{abstract}
\tableofcontents
\section{问题重述}
\subsection{问题的提出}
\section{模型的假设}
这里有一双小眼睛

\centerline{\includegraphics[width=1cm]{eye.png}}
\section{符号说明}
\begin{center}
\begin{tabular}{cc}
\hline
\makebox[0.3\textwidth][c]{符号} & \makebox[0.4\textwidth][c]{意义} \\ \hline
D & 木条宽度（cm） \\ \hline
\end{tabular}
\end{center}
\section{问题分析}
\section{总结}
\begin{thebibliography}{9}%宽度9
\bibitem{bib:one} ....
\end{thebibliography}
\begin{appendices}
附录的内容。
\end{appendices}
\end{document}

你就会得到：看着还是很清爽的

 
   当然“熟读唐诗三百首，不会作诗也能吟”。想要写得好论文，首先是要多看往年的好论文。
   

写在最后

   虽然目前仍算对建模“一无所知”，但相信在走向建模的过程，收获将会一点一滴累积。

   所以，好好加油，day day up!
   

---

1. 参照卓金武、王鸿钧的《MATLAB数学建模方法与实践(第3版)》中框架 ↩︎

2. 算法程序整理自网络，下同 ↩︎