建模方法笔记

本文简要介绍了一些数学建模常用知识点，以及相关的数学公式，仅供参考，每个知识点或许不太深入，后续有机会将会一个点一个点完善。

线性规划问题

定义：目标函数+决策变量+约束条件，目标函数与约束条件均为线性函数

一般线性规划问题的标准型：
$$
\max z=\sum_{j=1}^nc_jx_j\
s.t.\begin{cases}
\sum_{j=1}^n{a_{ij}x_j=b_i} \enspace i=1,2,···,m\
x_j\ge 0 \enspace j=1,2,···,n\
\end{cases}\
$$
可行解：满足一组约束条件的解

最优解：使目标函数达到最大值的可行解

可行域：所有可行解构成的集合，记为R

matlab中：
$$
\min\limits_{x} c^Tx\
s.t.\begin{cases}
Ax\leq b\
Aeq·x=beq\
lb\leq x \leq ub
\end{cases}\
其中c和x为n维向量，A、Aeq为适当维数的矩阵，b、beq为适当维数的列向量
$$
matlab求解：首先变为matlab标准型

[x,fval]=linprog(c,A,b)
[x,fval]=linprog(c,A,b,Aeq,beq)
[x,fval]=linprog(c,A,b,Aeq,beq,lb,ub)
% x返回的使决策向量的取值，fval返回的使目标函数的最优值，c为价值向量，A，b对应的使线性不等式约束，Aeq,beq对应的使线性等式约束，lb,ub分别对应的决策向量的下界和上界向量

可转化为线性规划的问题

$$
\min |x_1|+|x_2|+···+|x_n|\
s.t.\enspace Ax\le b\
对x_i存在任意的u_i,v_i>0满足\
x_i=u_i-v_i,|x_i|=u_i+v_i\
即取u_i=\frac{x_i+|x_i|}{2},v_i=\frac{|x_i|-x_i}{2}\
将原式变为：\min \sum_{i=1}^n(u_i+v_i)
\s.t.\begin{cases}
A(u-v)\leq b\
u,v\ge 0
\end{cases}
$$

整数规划模型

纯整数规划：所有决策变量要求取非负整数(引进的松弛变量和剩余变量可以不要求取整)

混合整数规划：只有一部分的决策变量要求取非负整数，另一部分可以取非负实数

全整数规划：除了所有决策变量要求取非负整数外，系数a_ij和常熟b_i也要求取整数

0-1整数规划：所有决策变量只能取0或1

特点：

原线性规划最优解全是整数，则整数规划最优解与线性规划最优解一致
整数规划可无可行解
有可行解，但最优解变差

分支定界法

先求出对应松弛问题的最优解
- 若松弛问题无可行解，则ILP无可行解
- 若松弛问题最优解符合整数要求，则就是ILP的最优解
- 若最优解不满足整数要求则构造新约束
选取不满足整数条件的变量 x_i⁰ 来构造，形成两个子问题
$$
x_i\le \lfloor x_i^o \rfloor\
x_i\ge \lfloor x_i^o \rfloor +1
$$
对两个子问题求解，观察是否有可行解
- 若有整数解，则存为拟定最优解，计算另一子问题的下两子问题是否存在更优解
- 若无重复第二步

分支定界法代码（matlab) _{文件Branch_Delimitation_methode}

% 主要修改参数即可进行分支定界的求解
clear;clc;
f = [-40, -90];
a = [8,7;7,20;];
b = [56;70];
aeq = [];
beq = [];
lb = [0; 0];
ub = [inf; inf];
Branch_Delimitation_method(f,a,b,aeq,beq,lb,ub)

% 分支定界法
function Branch_Delimitation_method(f, A, b, Aeq, beq, lbnd, ubnd)

global result; % 存储所有整数解
global lowerBound; % 下界
global upperBound; % 上界
global count; % 用以判断是否为第一次分支
count = 1; 

createBinTreeNode({f, A, b, Aeq, beq, lbnd, ubnd});
if ~isempty(result)
    [~,flag]=min(result(:,end)); % result中每一行对应一个整数解及对应的函数值
    Result=result(flag,:);
    disp('该整数规划问题的最优解为：');
    disp(Result(1,1:end-1));
    disp('该整数规划问题的最优值为：');
    disp(Result(1,end));
else
    disp('该整数规划问题无可行解');
end

BinTree函数

% 构建二叉树，每一结点对应一解 function BinTree = createBinTreeNode(binTreeNode) global result; global lowerBound; global upperBound; global count; if isempty(binTreeNode) return; else BinTree{1} = binTreeNode; BinTree{2} = []; BinTree{3} = []; [x, fval, exitflag] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, binTreeNode{7}); if count == 1 % upperBound = 0; % 初始下界为空 lowerBound = fval; count = 2; end if ~IsInRange(fval) return; end if exitflag == 1 flag = []; % 寻找非整数解分量 for i = 1 : length(x) if round(x(i)) ~= x(i) flag = i; break; end end % 分支 if ~isempty(flag) lowerBound = max([lowerBound; fval]); OutputLowerAndUpperBounds(); lbnd = binTreeNode{6}; ubnd = binTreeNode{7}; lbnd(flag, 1) = ceil(x(flag, 1)); % 朝正无穷四舍五入 ubnd(flag, 1) = floor(x(flag, 1)); % 进行比较，优先选择目标函数较大的进行分支 [~, fval1] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd); [~, fval2] = linprog(binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}); if fval1 < fval2 % 创建左子树 BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd}); % 创建右子树 BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}}); else % 创建右子树 BinTree{3} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, lbnd, binTreeNode{7}}); % 创建左子树 BinTree{2} = createBinTreeNode({binTreeNode{1}, binTreeNode{2}, binTreeNode{3}, ... binTreeNode{4}, binTreeNode{5}, binTreeNode{6}, ubnd}); end else upperBound = min([upperBound; fval]); OutputLowerAndUpperBounds(); result = [result; [x', fval]]; return; end else % 剪枝 return; end end

IsInRange剪枝函数和OutputLowerAndUpperBounds上下界输出函数

% 减枝 function y = IsInRange(fval) global lowerBound; global upperBound; if isempty(upperBound) & fval >= lowerBound y = 1; else if (fval >= lowerBound & fval <= upperBound) y = 1; else y = 0; end end % 打印输出上下界 function y = OutputLowerAndUpperBounds() global lowerBound; global upperBound; disp("此时下界、上界分别为"); disp(lowerBound); if isempty(upperBound) disp(' 无上界') else disp(upperBound); end end

割平面算法

先求出对应松弛问题的最优解
- 若松弛问题无可行解，则ILP无可行解
- 若松弛问题最优解符合整数要求，则就是ILP的最优解
若最优解不满足整数要求则增加割平面条件
$$
x_i+\sum a_{ik}x_k=b_i\enspace ,x_k 为松弛变量\
a_{ik}=[a_{ik}]+(a_{ik}-La_ik)=[a_{ik}]+f_{ik}\
b_i=[b_i]+(b_i-Lb_i)=[b_i]+f_i\
原式：x_i+\sum[a_{ik}]x_k+\sum f_{ik}x_i=[b_i]+f_i\
$$
增加一个线性约束，将松弛问题的可行区域割掉一块，使得非整数解恰好在割掉的一块中，但又没有割掉原问题的可行解，得到新问题
$$
移项，使得整数部分=小数部分\
x_i+\sum [a_{ik}]x_k -[b_i]=f_i-\sum f_{ik}x_i \leq 0\
f_i-\sum f_{ik}x_i\leq 0 \enspace 即为增加的线性约束
$$
若仍不为整数解则重复步骤

割平面法代码：_{文件Cut_plane_method}

function  [intx,intf] = DividePlane(A,c,b,baseVector)
%功能：用割平面法求解整数规划
%调用格式：[intx,intf]=DividePlane(A,c,b,baseVector)
%其中，A：约束矩阵；
%      c：目标函数系数向量；
%      b：约束右端向量；
%      baseVector：初始基向量；
%      intx：目标函数取最小值时的自变量值；
%      intf：目标函数的最小值；
sz = size(A);
nVia = sz(2);%获取有多少决策变量
n = sz(1);%获取有多少约束条件
xx = 1:nVia;

if length(baseVector) ~= n
    disp('基变量的个数要与约束矩阵的行数相等！');
    mx = NaN;
    mf = NaN;
    return;
end
 
M = 0;
sigma = -[transpose(c) zeros(1,(nVia-length(c)))];
xb = b;
 
%首先用单纯形法求出最优解
while 1   
    [maxs,ind] = max(sigma);
%--------------------用单纯形法求最优解--------------------------------------
    if maxs <= 0   %当检验数均小于0时，求得最优解。      
        vr = find(c~=0 ,1,'last');
        for l=1:vr
            ele = find(baseVector == l,1);
            if(isempty(ele))
                mx(l) = 0;
            else
                mx(l)=xb(ele);
            end
        end
        if max(abs(round(mx) - mx))<1.0e-7  %判断最优解是否为整数解，如果是整数解。
            intx = mx;
            intf = mx*c;
            return;
        else  %如果最优解不是整数解时，构建切割方程
            sz = size(A);
            sr = sz(1);
            sc = sz(2);
            [max_x, index_x] = max(abs(round(mx) - mx));
            [isB, num] = find(index_x == baseVector);
            fi = xb(num) - floor(xb(num));
            for i=1:(index_x-1)
                Atmp(1,i) = A(num,i) - floor(A(num,i));
            end
            for i=(index_x+1):sc
                Atmp(1,i) = A(num,i) - floor(A(num,i));
            end
            
            Atmp(1,index_x) = 0; %构建对偶单纯形法的初始表格
            A = [A zeros(sr,1);-Atmp(1,:) 1];
            xb = [xb;-fi];
            baseVector = [baseVector sc+1];
            sigma = [sigma 0];
         
            %-------------------对偶单纯形法的迭代过程----------------------
            while 1
                %----------------------------------------------------------
                if xb >= 0    %判断如果右端向量均大于0，求得最优解
                    if max(abs(round(xb) - xb))<1.0e-7   %如果用对偶单纯形法求得了整数解，则返回最优整数解
                        vr = find(c~=0 ,1,'last');
                        for l=1:vr
                            ele = find(baseVector == l,1);
                            if(isempty(ele))
                                mx_1(l) = 0;
                            else
                                mx_1(l)=xb(ele);
                            end
                        end
                        intx = mx_1;
                        intf = mx_1*c;
                        return;
                    else   %如果对偶单纯形法求得的最优解不是整数解，继续添加切割方程
                        sz = size(A);
                        sr = sz(1);
                        sc = sz(2);
                        [max_x, index_x] = max(abs(round(mx_1) - mx_1));
                        [isB, num] = find(index_x == baseVector);
                        fi = xb(num) - floor(xb(num));
                        for i=1:(index_x-1)
                            Atmp(1,i) = A(num,i) - floor(A(num,i));
                        end
                        for i=(index_x+1):sc
                            Atmp(1,i) = A(num,i) - floor(A(num,i));
                        end
                        Atmp(1,index_x) = 0;  %下一次对偶单纯形迭代的初始表格
                        A = [A zeros(sr,1);-Atmp(1,:) 1];
                        xb = [xb;-fi];
                        baseVector = [baseVector sc+1];
                        sigma = [sigma 0];
                        continue;
                    end
                else   %如果右端向量不全大于0，则进行对偶单纯形法的换基变量过程
                    minb_1 = inf;
                    chagB_1 = inf;
                    sA = size(A);
                    [br,idb] = min(xb);
                    for j=1:sA(2)
                        if A(idb,j)<0
                            bm = sigma(j)/A(idb,j);
                            if bm<minb_1
                                minb_1 = bm;
                                chagB_1 = j;
                            end
                        end
                    end
                    sigma = sigma -A(idb,:)*minb_1;
                    xb(idb) = xb(idb)/A(idb,chagB_1);
                    A(idb,:) = A(idb,:)/A(idb,chagB_1);
                    for i =1:sA(1)
                        if i ~= idb
                            xb(i) = xb(i)-A(i,chagB_1)*xb(idb);
                            A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:);
                        end
                    end
                    baseVector(idb) = chagB_1;
                end
              %------------------------------------------------------------
            end 
            %--------------------对偶单纯形法的迭代过程---------------------    
        end     
    else     %如果检验数有不小于0的，则进行单纯形算法的迭代过程
        minb = inf;
        chagB = inf;
        for j=1:n
            if A(j,ind)>0
                bz = xb(j)/A(j,ind);
                if bz<minb
                    minb = bz;
                    chagB = j;
                end
            end
        end
        sigma = sigma -A(chagB,:)*maxs/A(chagB,ind);
        xb(chagB) = xb(chagB)/A(chagB,ind);
        A(chagB,:) = A(chagB,:)/A(chagB,ind);
        for i =1:n
            if i ~= chagB
                xb(i) = xb(i)-A(i,ind)*xb(chagB);
                A(i,:) = A(i,:) - A(i,ind)*A(chagB,:);
            end
        end
        baseVector(chagB) = ind;
    end
    M = M + 1;
    if (M == 1000000)
        disp('找不到最优解！');
        mx = NaN; 
        minf = NaN;
        return;
    end
end

匈牙利算法

将系数矩阵转换为新矩阵
- 每行元素都减去改行的最小元素
- 每列元素都减去改列的最小元素
进行试指派，以寻求最优解
- 从只有一个零元素的行(列)开始，给他加圈
- 划去所在行(列)的其他零元素
- 给只有一个零元素的列(行)中的零元素加圈，划去加圈所在行的零元素
- 重复步骤
- 若仍有没有划圈的零元素，且同行（列）的零元素至少有两个，则从零元素最少的行（列）开始，比较所在列中的零元素个数，选择手的那个加圈
作最少的直线覆盖所有零元素
- 对没有圈的行打勾
- 对已打勾的行中所有含划去元素的列打勾
- 对有打勾的列中被圈出来的行打勾
- 重复步骤
对没有打勾的行画横线，有打勾的列画纵线，得到条数L
变换矩阵以增加零元素
- 从没有呗直线覆盖的所有元素中找到最小元素，打勾的每行每列都减去它

文件Hungarian_algorithm

非线性规划模型

一般形式：
$$
\min f(x)\
s.t.\begin{cases}
h_j(x)\le 0 \enspace j=1,2,····q\
g_i(x)=0 \enspace i=1,2,···p
\end{cases}
$$
matlab中：
$$
\min f(x)\
s.t.\begin{cases}
A·x \le b\
Aeq·x=beq\
c(x)\le 0\
ceq(x)=0\
lb\le x \le ub
\end{cases}\
c(x),ceq(x)是非线性向量函数
$$

[x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
fun:非线性函数，M文件定义的函数，目标函数
x0:x的初始值，随机在x中一个值
nonlcon:非线性约束条件，M文件定义的函数
options:优化参数，一般缺省[]

一般步骤：

编写函数fun1.m定义目标函数
- function f = fun1(x);
- f=sum(x.^2)+8;
编写函数fun2.m定义非线性约束条件
- function [g,h]=fun2(x);
- g=[] 非线性不等式约束
- h=[] 非线性等式约束
主程序文件

二次规划

定义：目标函数的自变量为二次函数，约束条件全为线性
$$
\min \frac{1}{2}x^THx+f^Tx\
s.t.\begin{cases}
Ax \le b\
Aeq·x=beq\
lb \le x \le ub
\end{cases}
$$

1
2
3

[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,option)
H:二次系数矩阵，平方处需要两倍
f:一次系数矩阵

微分方程

常微分方程

解析解
- 分离变量
- 一阶线性微分方程
数值解
$$
已知 \frac{dy}{dx}=e^{x^2},初始值：y(0)=1\
求 y = f(x)
$$
- 通过$\frac{f(x+h)-f(x)}{h}=e^{x^2}$ 得到递推公式$f(x+h)=f(x)+h\times e^{x^2}$
  1
  2
  3
  4
  5
  6
  7
  8
  clc;clear;
  k=1;x(1)=0;y(1)=1;
  for i=0:h:1
  y(k+1)=y(k)+exp(i^2)*h;
  x(k+1)=i+k;
  k=k+1;
  end
  plot(x,y)

高阶方程

$$
\frac{d^2y}{dx^2}+\frac{dy}{dx}+y=e^{x^2}
$$

换元 $z=\frac{dy}{dx}$ 即 $$
同时求解$\begin{cases}y(x+h)=y(x)+z(x)\times h\ z(x+h)=z(x)+[e^{x^2}-z(x)-y(x)]\times h \end{cases}$

clc;clear;
k=1;x(1)=0;y(1)=1;z(1)=1;
h=0.01;
for i=0:h:1
    z(k+1)=z(k)+(exp(i^2)-z(k)-y(k))*h;
    y(k+1)=y(k)+z(k)*h;
    x(k+1)=i+h;
    k=k+1;
end
plot(x,y)
hold on
plot(x,z)

Matlab求微分方程（组）解析解

1 2	dsolve('方程1','方程2',...,'方程n','初始条件','自变量') D表示微分，D2、D3表示求二阶、三阶

$\frac{dy}{dx}=1+y^2,y(0)=1$

输入：y1=dsolve(‘Dy=1+y^2’,’x’) 求通解

y2=dsolve(‘Dy=1+y^2’,’y(0)=1’,’x’) 求特解

Matlab求微分方程（组）数值解

[T,Y]=solver(odefun,tspan,y0)
y0为初始条件
tspan为求解区间
solver：ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb

时间序列

时间序列：按时间顺序排列的一组随机变量${X_t}$
观测值序列：时间序列的n个有序观测值${x_t}$

时间序列分析的主要目的：

揭示时间序列变化的统计规律性
预测未来事件
控制将来事件

特征统计量

均值：$\mu_X(t)=E[X_t], t\in T $ (期望)
方差：$D_X(t)=\gamma_X(t,t)=E[X_t-\mu_X(t)]^2,t\in T$
自协方差：$\gamma_X(t,s)=E[(X_t-\mu_X(t))(X_s-\mu_X(s))], t,s\in T$
自相关系数：$\rho_X(t,s)=\frac{\gamma_X(t,s)}{\sqrt{D_X(t)·D_X(s)}}$

特征统计量的估计

样本均值 $\hat{\mu}=\overline{x}=\frac{1}{n}\sum_{i=1}^n x_i$
样本方差 $\hat{\gamma}(0)=\frac{1}{n}\sum_{i=1}^n(x_i-\overline{x})^2$
样本延迟k自协方差 $\hat{\gamma}(k)=\frac{1}{n}\sum_{i=1}^{n-k}(x_i-\overline{x})(x_{i-k}-\overline{x})$
样本延迟k自相关系数 $\hat{\rho}(k)=\frac{\hat{\gamma}(0)}{\hat{\gamma}(k)}$

平稳性

严平稳：事件序列的概率分布与事件无关
宽平稳：$E[X_t],E[X_t^2]$为常数，$\gamma_X(t,s)=\gamma_X(k+t,k+s)$

模型分类

线性事件序列模型
- 自回归滑动平均（ARMA）模型
- 自回归综合滑动平均（ARIMA）模型
- 季节ARIMA（SARIMA）模型
非线性事件序列模型
- 自激励门限自回归（SETAR）模型
- 双线性（BL）模型
- 指数自回归（EAR）模型

ARMA模型

AR自回归模型
MA滑动平均模型
ARMA模型

AR(p) 自回归模型

一般形式：
$$
X_t=\phi_0+\phi_1x_{t-1}+…+\phi_px_{t-p}+\xi_t
$$
对于前期数据的影响

MA 滑动平均模型

一般形式
$$
X_t=\mu+\xi_t-\theta_1\xi_{t-1}-…\theta_q\xi_{t-q}
$$
不考虑历史数据，考虑扰动

ARMA模型

结合AR,MA
$$
X_t=\phi_0+\phi_1x_{t-1}+…+\phi_px_{t-p}+\xi_t+\mu+\xi_t-\theta_1\xi_{t-1}-…\theta_q\xi_{t-q}
$$

最小二乘法

给定参数
- 数据确定范围
  - 确定最大斜率和最小斜率
- 确定初始值
- 粗/细搜索法
计算误差
寻找最小误差

生长模型

对数生长模型：假设肿瘤增长速度与其大小的对数成正比，即每个时间单位内肿瘤的增长量与其当前大小的对数成比例。这种模型适合于描述肿瘤在较短时间内的生长情况。
Gompertz生长模型：该模型假设肿瘤生长速度与其大小的指数成负相关，即随着时间的推移，肿瘤的生长速度逐渐减缓。这种模型适用于描述肿瘤在较长时间内的生长情况。
Logistic生长模型：该模型假设肿瘤生长速度与其大小的二次函数成正相关，但是随着时间的推移，增长速度逐渐减缓，直到趋于稳定。这种模型适用于描述肿瘤在中期生长情况。
多阶段生长模型：该模型将肿瘤的生长过程分为多个阶段，并在每个阶段内采用不同的生长模型进行描述。这种模型适用于描述肿瘤生长过程中转变的情况。
机器学习模型：除了基于数学模型的建模方法外，还可以采用机器学习模型，如神经网络、支持向量机等方法，通过对大量实验数据的学习，建立肿瘤生长模型，并对肿瘤的生长情况进行预测。

Gompertz生长模型

$$
\frac{dx}{dt}=rx\ln\frac{k}{x}
$$

x为种群的数量
k为环境最大容纳量
r为增长率的系数

Logistic生长模型

$$
\frac{dx}{dt}=rx
$$

x为种群数量
r为增长率系数

阻滞增长

用r(x)表示r关于x的变化:r(x)=r-sr
$$
\frac{dx}{dt}=r(1-\frac{x}{x_m})x
$$
解析解为：
$$
x(t)=\frac{k}{1+(\frac{k}{x_0}-1)e^{-rt}}
$$

放疗系数

以及放疗对于肿瘤细胞的增长有着抑制作用,我们用$\alpha$来描述放疗的杀伤效应，可以得到模型：
$$
\begin{cases}
\frac{\text{d}x}{\text{d}t}=r(1-\frac{x}{k})x-\alpha x\
\
x(0)=x_0
\end{cases}
$$
假设放疗的杀伤作用是比较稳定的，即放疗对细胞的杀伤率与时间成正比（假设5.），我们使用指数函数来描述它$\alpha = \alpha_0 e ^{-\beta t}$

$$
\begin{cases}
\frac{\text{d}x}{\text{d}t}=r(1-\frac{x}{k})x-\alpha_0e^{-\beta t} x\
\
x(0)=x_0
\end{cases}
$$
解析解为：
$$
x(t)=\frac{k}{(1+\alpha e^{-\beta t})(1+(\frac{k}{x_0}-1)e^{-rt})}
$$

综合评价

构成综合评价问题的五要素

被评价对象、评价指标、权重系数、综合评价模型、评价者

评价指标：系统性、科学性、可比性、可测性、独立性

综合评价的一般步骤

明确评价目的，确定被评价对象
建立评价指标系（评价指标的原始值、评价指标的若干预处理）
确定与各项评价指标相对应的权重系数
选择或构造综合评价模型
计算系统的综合评价值并给出评价结果

指标规范化

指标类型正向化
标准化处理
- 标准差
- 极值差法
- 功效系数法

赋权法

基于指标差异的赋权法
- 均方差法
- 极差法
- 熵值法
动态赋值方法

层次分析法基本模型（AHP）

特点：对复杂的决策问题的本质、影响因素以及其内在关系等进行深入分析的基础上，利用较少的定量信息使决策的思维过程数学化，从而为多目标、多准测或无结构特性的复杂决策问题提供简便的决策方法

**文件 **Analytic_hierarchy

三大典型应用

用于最佳方案的选取
用于评价类问题
用于指标体系的优化

基本原理

将问题分解为不同的组成因素，并按照因素间的相互关联影响以及隶属关系，将因素按不同层次聚集组合，形成一个多层次的分析结构模型，将问题归结为最低层相对于最高层的相对重要权值的确定或相对优劣次序的排定

建立层次结构模型
构造判断（成对比较）矩阵
层次单排序以其一致性检验
层次总排序及其一致性检验

基本过程

建立层次结构模型
- 最高层：目标层，决策的目的、要解决的问题
- 中间层：准测层，考虑的因素、决策的准则
- 最底层：方案层，决策时的备选方案
- 对于相邻的两层，高层为目标层，低层为因素层
成对比较
- 不把所有因素放在一起比较，而是两两比较
- 采用相对尺度，减少性质不同的诸多因素相互比较的困难
- 使用1~9的标度方法
  
  标度含义
  
  1 同样重要
  
  3 稍微重要
  
  5 明显重要
  
  7 强烈重要
  
  9 极端重要
  
  2，4，6，8 上述两相邻判断的中值
  
  倒数不重要程度
- 构建判断矩阵（一致矩阵or不一致矩阵）
层次单排序及其一致性检验
- 求最大特征根$\lambda_{max}$ 以及特征向量$Aw=\lambda w$，经过归一化后记为 $w$，且$w$ 的元素为同一层次因素对于上一层次因素某因素相对重要性的排序权值
- 运用定理检验一致性
  - n阶一致阵的唯一非零特征根为n
  - n阶正互反阵A的最大特征根$\lambda \ge n$,当且仅到$\lambda = n$时A为一致阵
  - 定义一致性指标：$CI = \frac{\lambda-n}{n-1}$ 越接近0，越一致
- 计算CR指标：$CR = \frac{CI}{RI}$,$CR < 0.1$则在范围内，RI为随机一致性指标，查表可得
- 若为一致阵，则列向量归一化，求行和归一化，即可得到$w$
层次总排序及其一致性检验
- 构建每一指标对于方案的比较矩阵（几个指标几个矩阵）
- 方案得分：矩阵乘上层的权重求和
- $CR = \frac{CI_1+CI_2+···+CI_n}{RI_1+RI_2+···+RI_n} <0.1$

标度	含义
1	同样重要
3	稍微重要
5	明显重要
7	强烈重要
9	极端重要
2，4，6，8	上述两相邻判断的中值
倒数	不重要程度

局限性

粗略、主观，比较、判断以及结果都是粗糙的不适于精度要求很高的问题。其次是从建立层次结构图到给出两两比较矩阵，人的主观因素作用很大，使决策结果较大程度的依赖于决策人的主观意志，可能难以为众人所接收

TOPSIS模型算法

多指标评价方法，通过构造评价问题的正理想解和负理想解，即各指标的最优解和最劣解，通过计算每个方案到理想方案的相对贴近度，排序选择最优方案

文件： Evaluation –>TOPSIS(手动加权)

方法和原理

设多属性决策方案集 $D = {d_1,d_2,···,d_m}$ ,衡量方案优劣的属性变量为 $x_1,x_2,···,x_n$ ,这时方案集D中的每个方案$d_i$ 的n个属性值构成的向量是$[a_{i1},a_{i2},···,a_{in}]$ 。

正理想解，方案集D中并不存在的虚拟的最佳方案；负理想解，方案集D中并不存在的虚拟的最劣方案。

算法步骤

数据预处理（正向化，非量纲化，归一化）

线性变化
$$
b_{ij}=a_{ij}/a_j^{max} ,效益型\
b_{ij}=1-a_{ij}/a_j^{max},成本型
$$
标准0，1变化
$$
b_{ij}=\frac{a_{ij}-a_j^{min}}{a_j^{max}-a_j^{min}}，效益型\
b_{ij}=\frac{a_j^{max}-a_{ij}}{a_j^{max}-a_j^{min}}，成本型
$$

区间型属性变化

$$
b_{ij}=\begin{cases}
1-\frac{(a_j^0-a_{ij})}{a_j^0-a_j^{down}}, a_j^{down} \le a_{ij} \le a_j^0\
1, a_j^0 \le a_{ij} \le a_j^{up}\
1-\frac{a_{ij}-a_j^*}{a_j^{up}-a_j^*},a_j^* \le a_{ij} \le a_j^{up}\
0,else\enspace cases
\end{cases}\
a_j^{down}可接收的下线，一般取min\
a_j^{up}可接受的上线
$$

x_new = @(interval,lb,ub,x)(1-(interval(1)-x)./(interval(1)-lb)).*(x>=lb&x<interval(1))+(x>=interval(1)&x<=interval(2))+(1-(x-interval(2))./(ub-interval(2))).*(x>interval(2)&x<=ub);
interval=[];lb= ;ub= ;
x_data=[];
x=x_new(interval,lb,ub,x_data)

向量归一化（最常用,但不适合区间型）
$$
效益型b_{ij}=\frac{a_{ij}}{\sqrt{\sum_{i=1}^m a_{ij}^2}} \enspace ,i=1,2···，m,j=1,2···,n\
成本型b_{ij}=1-\frac{a_{ij}}{\sqrt{\sum_{i=1}^m a_{ij}^2}} \enspace ,i=1,2···，m,j=1,2···,n
$$
标准化处理zscore(x)
$$
\overline{a_j}=\frac{1}{m}\sum_{i=1}^ma_{ij}\
S_j = \sqrt{\frac{1}{m-1}\sum_{i=1}^m(a_ij-\overline{x_j})^2}\
b_{ij}=\frac{a_{ij}-\overline{a_j}}{S_j}
$$

用向量规划的方法球的规范决策矩阵
- 消除纲量，使用决策矩阵$A = (a_{ij}){m\times n} $ ,正向化，规范化后的决策矩阵$B=(b{ij})_{m\times n} $
构成加权规范矩阵$C= (c_{ij})_{m\times n} $
- 赋予各属性权重 $w = [w_1,w_2,···,w_n]^T$
- $c_{ij}=w_j · b_{ij},i=1,2,···,m,j=1,2,···,n$
确定正理想解和负理想解
- 正理想解 $c_{j}^*$
- 负理想解 $c_{j}^0$
计算距离
- 计算备选方案$d_i$ 到正负理想解的距离
  $$
  S_i^* = \sqrt{\sum_{j=1}^n(c_{ij}-c_j^*)^2},i=1,2,···,m\
  S_i^0 = \sqrt{\sum_{j=1}^n(c_{ij}-c_j^0)^2},i=1,2,···,m\
  $$
计算最终指标
- $f_i^*=s_i^0/s_i^0+s_i^*) ,i=1,2···m$

熵权法

利用信息熵计算各个指标的权重，为多指标综合评价提供依据

归一化（消除量纲的影响）
$$
\frac{x-min}{max-min}(min-max归一化)
$$
计算信息熵
$$
H(p)=-\sum_{i=1}^np_i\log_2p_i
$$
p为概率：$\frac{元素的值}{列和}$

计算每一列的信息熵：
$$
e_j=-\frac{1}{\ln n}\sum_{i=1}^n[p_{ij}\times \ln(p_{ij})]\enspace,j=1,2,3…m
$$
确定权重
$$
d_j=1-e_j\w_j=\frac{d_j}{\sum_{j=1}^md_j}\enspace,j=1,2,3…m
$$

聚类分析

将一个数据集划分为若干组或类的过程，并使得同一个组内的数据对象具有较高的相似度，而不同组中的数据对象是不相似的，通常利用各数据对象间的距离来表示。适用于探讨样本之间的互相关联关系从而对一个样本结构做一个初步评价

距离

欧氏距离$d(x_i,x_j)=[\sum_{k=1}^p(x_{ik}-x_{jk})^2]^\frac{1}{2}$ pdist(x)
绝对距离$d(x_i,x_j)=\sum_{k=1}^p|x_{ik}-x_{jk}|$ pdist(x,’cityblock’) 一般不使用
明氏距离$d(x_i,x_j)=[\sum_{k=1}^p|x_{ik}-x_{jk}|^m]^\frac{1}{m}$ pdist(x,’minkowski’,r)
切氏距离$d(x_i,x_j)=\max\limits_{1 \le k \le p}|x_{ik}-x_{jk}|$ max(abs(xi-xj))
方差加权距离$d(x_i,x_j)=[\sum_{k=1}^p(x_{ik}-x_{jk})^2/s_k^2]^\frac{1}{2}$ 讲原数据标准化以后的欧式距离
马氏距离$d(x_i,x_j)=\sqrt{(x_i-x_j)^T\sum^{-1}(x_i-x_j)}$ pdist(x,’mahal’)
兰氏距离$d(x_i,x_j)=\frac{1}{p}\sum_{k=1}^p\frac{|x_{ik}-x_{jk}|}{x_{ik}+x_{jk}}$
杰氏距离$d(x_i,x_j)=[\sum_{k=1}^p(\sqrt{x_{ik}}-\sqrt{x_{jk}})^2]^\frac{1}{2}$

距离矩阵

a=[];
d1=pdist(a); %欧氏距离举例
D=squareform(d1); %实对称矩阵
S=tril(D); %保留下三角

相似系数

当对p个指标变量进行聚类时，用相似系数来衡量变量之间的相似程度 $C_{\alpha.\beta}$ ，应该满足
$$
|C_{\alpha,\beta}|\le 1 ,且C_{\alpha.\alpha}=1\
C_{\alpha,\beta}= \pm1,当且仅当\alpha=k\beta\
C_{\alpha,\beta}=C_{\beta,\alpha}
$$
夹角余弦
$$
C_{ij}(1)=\cos \alpha_{ij}=\frac{\sum_{t=1}^n x_{ti}x_{tj}}{\sqrt{\sum_{t=1}^nx_{ti}^2}\sqrt{\sum_{t=1}^nx_{tj}^2}}
$$
相关系数
$$
C_{ij}(2)=\frac{\sum_{t=1}^n(x_{ti}-\overline{x_i})(x_{tj}-\overline{x_j})}{\sqrt{\sum_{t=1}^n(x_{ti}-\overline{x_i})^2}\sqrt{\sum_{t=1}^n(x_{tj}-\overline{x_j})^2}}
$$

1
2
3

R=corrcoef(a); %指标之间的相关系数
a1=normc(a); %将a的各列化为单位向量
J=a1'*a1 %计算各列之间的夹角余弦

类间距离

$d_{ij}$ 表示两个样本 $x_i,x_j$ 之间的距离，$G_p,G_q$ 分别表示两个类别，各自含有 $n_p,n_q$ 个样本

最短距离 $D_{pq} = \min\limits_{i\in G_p,j \in G_q} d_{ij}$ linkage(d)
最长距离 $D_{pq} = \max\limits_{i\in G_p,j \in G_q} d_{ij}$ linkage(d,’conplete’)
重心距离 $D_{pq} = d(\overline{x_p},\overline{x_q})= \sqrt{(\overline{x_p}-\overline{x_q})^T(\overline{x_p}-\overline{x_q})}$ linkage(d,’average’)
类平均距离 $D_{pq} = \frac{1}{n_pn_q}\sum_{i\in G_p}\sum_{j \in G_q}d_{ij}$ linkage(d,’centroid’)
离差平方和距离 $D_{pq}^2=\frac{n_pn_q}{n_p+n_q}(\overline{x_p}-\overline{x_q})^T(\overline{x_p}-\overline{x_q})$ linkage(d,’ward’)

谱系聚类法

选择样本间距离的定义以及类间距离的定义
构造个类，每类只含有一个样本,计算n个样本两两间的距离，得到距离矩阵D
- n个样本开始作为n个类，计算两两间的距离或是相似系数，得到实对称矩阵$D_0$
合并符合个类间距离定义要求的两类为一个新类
- 从$D_0$的非主对角线上找最小或是最大的元素，设为$D_{pq}$ ,将$G_p,G_q$ 合并为一个新类，并从$D_0$中去除两行两列，在加上新类与个类之间的距离，得到一个n-1阶的矩阵
计算新类与当前各类的距离，若类的个数为1则第六步，否则第四步
画出聚类图
- H=dendrogram(z,d) %样本少于30可以不写d
- T = cluster(z,k) % k为分类数目
- Find(T==k0) % 找出属于第k0类的样品编号
决定类的个数和类

K-平均聚类算法（K-means）

从n个数据对象任意选择k个对象作为初始聚类中心
根据每个聚类对象的均值，计算每个对象到这些中心对象的距离，并根据最小距离重新划分
重新计算每个类的均值
循环直至不在变化

特点： 只适用于聚类均值有意义的场合，对噪声和孤立点数据敏感

灰色预测模型

数据少，看不出明显规律，适合灰色预测

制造规律
- 累加生成序列$x^{(1)}(k)=\sum_{i=1}^kx^{(0)}(i)$
- 再修正为均值生成序列，使用前后两个时刻的均值
新序列像指数曲线
- 使用一个指数曲线乃至一条直线逼近序列
- 构建一阶常微分方程求解拟合曲线
设 $x^{(1)}$ 满足$\frac{dx^{(1)}}{dt}+ax^{(1)}=u$
- 使用最小二乘法求a和u
数据是离散的
- $\frac{dx^{(1)}}{dt}=\frac{\Delta x^{(1)}}{\Delta t}$
- $\Delta t =(t+1)-t=1$ 始终为1
- $\Delta x^{(1)}=x^{(1)}(t)-x^{(1)}(t-1)=x^{(0)}(t)$
- 由此得到 $x^{(0)}(t)=-ax^{(1)}(t)+u$
模型检验
- 残差检验：（真实-预测）/ 预测，以0.2为界限

数据预处理

数据清洗
- 去除掉数据中的噪声，纠正不一致
- 缺失值处理：删除记录、数据插补、不处理
- 插补方法：
  - 均值/中位数/众数
  - 固定值
  - 最近临插补法
  - 回归
  - 插值法
- 插值法：
  - 拉格朗日插值法（与插值点个数有关
  - 牛顿插值法（与插值点个数无关
数据集成
- 将多个数据源合并成一致的数据存储，构成一个完整的数据集
数据变换
- 归一化
- 最大最小规范化：$x^*=\frac{x-min}{max-min}$
- 零均值规范化：$x^*=\frac{x-\overline{x}}{\sigma},其中\overline{x}为均值，\sigma为标准差$
数据规约
- 通过集成、删除冗余属性或聚类等方法来压缩数据

判断是否为异常值

计算n为数据集中所有样本间的测量距离
若样本S中之手有一部分数量为P的样本到$S_i$的距离比d大，则为噪声数据

特征工程

特征选取
- 过滤法
- 包装法
- 嵌入法
特征预处理

插值和拟合算法

插值

拉格朗日插值 lagrange(x,y,x0)

$f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)$ 每一次对一个点拟合，最后相加

function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
    z=x(i);
    s=0.0;
    for k=1:n
        p=1.0;
        for j=1:n
            if j~=k
                p=p*(z-x0(j))/(x0(k)-x0(j));
            end
        end
        s=p*y0(k)+s;
    end
    y(i)=s;
end
1
2
3
4
5
6
7
8
9
10
11
12

    

- 高次插值的Runge现象

  - 插值多项式的次数超过七次，会出严重的振荡现象
  - 避免Runge现象的常用方法：分段插值

#### matlab插值算法

- 一维插值interp1


yi=interp1(x,y,xi,'method')
nearest:最近邻插值
linear:线性插值
spline:三次样条插值,一般优先
cubic:立方插值
省略为线性插值
1


x=0:2:24;
y-[12 9 9 10 18 24 28 27 25 20 18 15 13];
x1=13;
y1=interp1(x,y,x1,'spline')
xi=0:1/3600:24;
yi=interp1(x,y,xi,'spline');
plot(x,y,'*',xi,yi)
1
2
3

- 二维插值interp2


zi=inter2(x,y,z.xi,yi,'method')
nearest:最邻近插值
linear:双线性插值
spline:双三次样插值
cubic:双立方插值，一般优先
1


x=1:5;
y=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
figure(1);
mesh(x,y,temps);
xi=1:0.2:5;
yi=1:0.2:3;
zi=interp2(x,y,temps,xi.yi,'cubic');
figure(2);
mesh(xi,yi,zi);
figure(3);
contour(xi,yi,zi,20,'r');	% 等高线
[i,j]=find(zi==min(min(zi)));
x=xi(j),y=yi(i),zmin=zi(i,j)
[i,j]=find(zi==max(max(zi)));
x=xi(j),y=yi(i),zmax=zi(i,j)
1
2
3
4
5

  plot3空间曲线，mesh空间曲面，surf空间曲面，contour等高线

- 散乱点插值griddata


griddata(x,y,z,xi,yi,'method')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15


插值主要用于求函数值，拟合主要求函数关系

#### 拟合

1. 线型的选择
   - 若是线性，采用最小二乘法
   - 若是非线性，采用Gauss-Newton迭代法
2. 线型中参数的计算

#### Matlab拟合

- 多项式拟合


[a,S]=polyfit(x,y,n)
n为拟合多项式的次数
a为多项式系数构成的向量
S为分析拟合效果所需的指标
1


x=1:12;
y=[5 8 9 15 25 29 31 30 22 25 27 24];
a=polyfit(x,y,9)
xp=1:0.1:12;
yp=polyval(a,xp);
plot(x,y,'.k',xp,yp,'r');
1
2
3

- 非线性拟合


[b,r]=polyfit(x,y,fun,b0,option)
fun为拟合函数
b0为拟合参数的初始迭代值
option为拟合选项
b为拟合参数
r为拟合残差
1
2
3
4
5
6
7
8
9
10

```;
x=1:16;
y=[4.00 6.40 8.00 8.80 9.22 9.50 9.70 9.86 10.00 10.20 10.32 10.42 10.50 10.55 10.58 10.60];
y1=@(b,t)b(1)*exp(-t/b(2))+b(3)*exp(-t/b(4))+b(5);
b0=[-1 1 -1 1 1];
a=nlinfit(x,y,y1,b0)
xp=1:0.1:16;
yp=y1(a,xp);
plot(x,y,'.k',xp,yp,'r');

cftool工具箱

回归模型

一元线性回归

一元线性回归模型
$$
\begin{cases}
y=\beta_0+\beta_1x+\xi\
E_\xi=0,D_\xi=\sigma^2\
\end{cases}
$$
$\beta_0、\beta_1$为回归系数，自变量$x$为回归变量，$Y=\beta_0+\beta_1x$为$y$对$x$的回归直线方程
主要任务：
1. 用试验值（样本值）对$\beta_0,\beta_1,\sigma$作点估计
2. 对回归系数做假设检验
3. 在$x=x_0$处对$y$做预测，对$y$做区间估计
检验：
- 对回归方程$Y=\beta_0+\beta_1x$显著性检验，归结为对$H_0:\beta_1=0;H_1:\beta \neq0$ 进行检验，一般小于百分之五
- F检验法（假设性检验）
  - 当$H_0$成立时，$F=\frac{U}{Q_e/(n-2)}$~$F(1,n-2)$ 其中$U=\sum_{i=1}^n(\hat{y_i}-\overline{y})^2$ （回归平方和）
  - $F>F_{1-\alpha}(1,n-2)$,拒绝$H_0$ (查表)，拒绝说明效果好
- t检验法（<30 适用于样本含量较小）
  - 当$H_0$成立时，$T=\frac{\sqrt{L_{xx}}\hat{\beta_1}}{\hat{\sigma_e}}$~$t(n-2)$ 其中$L_{xx}=\sum_{i=1}^n(x_i-\overline{x})^2=\sum_{i=1}^nx^2-n\overline{x}^2$
  - $|T|>t_{1-\frac{\alpha}{2}}(n-2)$,拒绝$H_0$
  - $\sigma_e^2=Q_e/(n-2)$ $Q_e=\sum_{i=1}^n(y_i-\hat{y_i})^2$
- r检验法（相关系数检验）
置信区间
预测与控制

蚁群算法

图论问题

参数初始化
- 设蚂蚁数量为m，城市数量n，城市i和城市j之间的距离为d_ij
- t 时刻城市 i 和城市 j 之间的路径上信息素浓度$\tau_{ij}(t)$
- 初始时各路径的信息素浓度相同，设为$\tau_{ij}(0)=\tau_0$
计算概率
- 当前蚂蚁再城市i ,假设只考虑信息素浓度对蚂蚁选择路径的影响：
  $$
  P_{ij}^k(t)=\begin{cases}
  \frac{[\tau_{ij}(t)]^\alpha \times [\frac{1}{d_{ij}(t)}]^\beta}{\sum_{s\in allowed_k}[\tau_{is}(t)]^\alpha \times [\frac{1}{d_{is}(t)}]^\beta },j\in allowed_k\
  0, j\notin allowed_k
  \end{cases}\
  allowed_k:第k只蚂蚁暂未访问的城市集合\
  P_{ij}^k:t时刻第k只蚂蚁从i到j的概率\
  s:暂未访问的城市的集合中的某一城市\
  p_{ij}^k的值越大，前往j的概率越大，但是代码中使用
  $$
更新信息素浓度
$$
\begin{cases}
\tau_{ij}(t+1)=(1-\rho)\tau_{ij}(t)+\sum_{k=1}^m \Delta \tau_{ij}^k, 0<\rho<1\
\Delta \tau_{ij}^k = \begin{cases}
\frac{Q}{L_k},第k只蚂蚁曾经过路径i到j\
0,其他
\end{cases}
\end{cases}\
$$
判断是否终止
- 如果达到最大迭代次数，算法终止，输出最优解

贪心算法

采用贪心思想，保证每次操作都是局部最优的，从而使最后得到的结果是全局最优的

弗洛伊德算法

图论中求解最短路径

蒙特卡洛

随机模拟方法，使用随机数解决问题，得到的解都是近似解

蒙特卡洛的基本原理简单描述是先大量模拟，然后计算一个事件发生的次数，再通过这个发生次数除以总模拟次数，得到想要的结果。比如投3个骰子，计算3个骰子同时是6的概率，可以模拟投N次（随机样本数），统计同时是6出现的次数C，然后C除以N即是计算结果。