数学建模
1.层次分析法
解决问题
评价类问题，即选哪种方案更好
基本原理
考虑影响结果的不同指标（这个指标可以通过知网查文献获取），并且知道它们的重要程度，构造一个判断矩阵
在用判断矩阵求权重之前，先进行一致性检验
引理：\(n\) 阶正反矩阵 \(A\)为一致矩阵时，当且仅当最大特征值\(\lambda(max)=n\)
为非一致矩阵时，\(\lambda(max)>n\)
越不一致，与 \(n\) 相差越大
一致性检验步骤：
1.计算一致性指标
2.查找对应的平均随机一致性指标
3.计算一致性比例
同时我们还需要计算权重，有三种方法，检验通过权重时才能够继续用
代码
1 构建判断矩阵
点击查看代码
import pandas as pd
import numpy as np
from functools import reduce
#定义两个元素相乘
import math
def chf(x,y):
return x*y
###1.构建判断矩阵
#1.1方法一，手动录入全部矩阵
matrix_initial= np.array([[1,0.333,1,3],[3,1,1,3],[1,1,1,3],[0.333,0.333,0.333,1]])
print(matrix_initial)
#1.1.1取得行数和列数
[m,n] = a.shape
或者
点击查看代码
import pandas as pd
import numpy as np
from functools import reduce
#定义两个元素相乘
import math
def chf(x,y):
return x*y
#1.2方法二，部分录入矩阵内容
#1.2.1构建对角线为1的矩阵
n = int(input("请输入维度数n，代码将根据输入数据构建n维矩阵,按回车键完成录入"))
matrix_initial = np.eye(n)
#1.2.2填写矩阵值
for j in range(0,n-1):
for i in range(j,n-1):
i=i+1
matrix_initial[i,j]=float(input("请输入a%s%s处矩阵值,按回车键完成录入"%(i,j)))
matrix_initial[j,i]=round(1/matrix_initial[i,j],6)
if matrix_initial[j,i]>1:
matrix_initial[j,i]=int(matrix_initial[j,i])
print("当前结构为\n%s"%matrix_initial)
j=j+1
print("判断矩阵最终为\n%s"%(matrix_initial))
2 方根法计算权重
点击查看代码
#2.方根法计算权重
#2.1每一行元素相乘并求出对应立方根，并存入列表。f为对应立方根,
matrix_mid = []
for i in range(0,n):
f = reduce(chf,matrix_initial[i,:])
f = pow(f,1/n)
matrix_mid.append(f)
matrix_mid = np.array(matrix_mid)
matrix_mid = matrix_mid.reshape(-1,1)
#计算立方根之和
matrix_weight = []
matrix_add = sum(matrix_mid)
#计算权重
for i in range(0,n):
f = matrix_mid[i,:]/matrix_add
matrix_weight.append(f)
matrix_weight = np.array(matrix_weight)
print(matrix_weight)
3 一致性检验
点击查看代码
#3.一致性检验
#特征值和特征向量
E,F = np.linalg.eig(matrix_initial)
#最大特征值
eigenval = np.max(E)
#CI和RI数据获取
CI=(eigenval-n)/(n-1)
RI=[0,0,0.58,0.9,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-1]
if CR>=0.1:
print('没有通过一致性检验\n')
else:
print('通过一致性检验\n')
print(CR)
print(matrix_weight)
参考网址
https://zhuanlan.zhihu.com/p/181711354
2.元胞自动机
解决问题
流行病的传播，生命游戏等
基本原理
每一个元胞下一时刻的状态只受到此时元胞自身和周围元胞的影响，全局变化是由一个个个体变化引起的
边界条件：
定值型：边界均为某个定值
周期型：像轮胎一样首尾相接
反射型：到边界外时呈现镜面反射
吸收型（绝热型）：边界外元胞和边界元胞的状态保持一致
代码
模拟森林火灾
点击查看代码
close;
clear;
clc;
n = 300; %元胞矩阵大小
Plight = 0.000001; Pgrowth = 0.001;
UL = [n 1:n-1];
DR = [2:n 1];
veg = zeros(n,n); %初始化
% The value of veg:
% empty == 0
% burning == 1
% green == 2
imh = image(cat(3,veg,veg,veg));
m=annotation('textbox',[0.1,0.1,0.1,0.1],'LineStyle','-','LineWidth',1,'String','123');
for i = 1:100000
sum = (veg(UL,:) == 1) + (veg(:,UL) == 1) + (veg(DR,:) == 1) + (veg(:,DR) == 1);
%根据规则更新森林矩阵：树 = 树 - 着火的树 + 新生的树
veg = 2 * (veg == 2) - ( (veg == 2) & (sum > 0 | (rand(n,n) < Plight)) ) + 2 * ( (veg == 0) & rand(n,n) < Pgrowth);
a=find(veg==2);
b=find(veg==1);
aa=length(a);
bb=length(b);
shu(i)=aa;
fire(i)=bb*30;
if (bb>=0&&bb<=10)
str1='森林正常';
elseif (bb>10&&bb<=100)
str1='火灾发展';
elseif (bb>100)
str1='森林大火';
end
if ((aa>48000)||(bb>=10))
str2='火灾预警：红色预警';
elseif (aa>42000&&aa<=48000)
str2='火灾预警：黄色预警';
elseif (aa>35000&&aa<=42000)
str2='火灾预警：蓝色预警';
elseif (aa>=0&&aa<=35000)
str2='火灾预警：安全';
end
str=[str1 10 str2];
set(imh, 'cdata', cat(3, (veg == 1), (veg == 2), zeros(n)) )
drawnow
figure(2)
delete(m)
plot(shu);
hold on
plot(fire);
legend(['绿树的数量',num2str(aa)],['火的数量',num2str(bb)]);
title(['时间T=',num2str(i),'天']);
m=annotation('textbox',[0.15,0.8,0.1,0.1],'LineStyle','-','LineWidth',1,'String',str);
hold off
% pause(0.0001)
end
参考网址
https://zhuanlan.zhihu.com/p/602124985
https://blog.csdn.net/qq_42112448/article/details/109514612
3 模拟退火
解决问题
在一个大的范围内寻找最优解（避免了像爬山算法一样陷入局部最优解，但也只是一定概率找到最优解）
基本原理
在初始阶段，尽可能地尝试不同的方向，到后期逐渐缩小，最终趋于稳定
初始化：
1.随机选择一个初始点 \(x_{0}\)
2.设置初始 “温度”，即控制参数 \(T_{0}\)
3.确定终止温度 \(T_{f}\) 达到这个温度或以下时，算法即可终止
如果 \(x_{new}\) 比 \(x\) 更好，就接受 \(x_{new}\)，如果 \(x_{new}\) 比不上 \(x\)，我们也有一定的概率接受 \(x_{new}\)，这个概率受到时间的影响
代码
点击查看代码
作者：小树
链接：https://zhuanlan.zhihu.com/p/676926594
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。
from simanneal import Annealer
import random
import numpy as np
import matplotlib.pyplot as plt
class TravelingSalesmanProblem(Annealer):
def move(self):
a = random.randint(0, len(self.state) - 1)
b = random.randint(0, len(self.state) - 1)
self.state[a], self.state[b] = self.state[b], self.state[a]
def energy(self):
e = sum(
dist_matrix[self.state[i-1]][self.state[i]]
for i in range(city_num)
)
return e
# 初始解（即初始城市访问顺序），随机产生
init_state = list(range(city_num))
random.shuffle(init_state)
tsp = TravelingSalesmanProblem(init_state)
tsp.set_schedule(tsp.auto(minutes=.2))
state, e = tsp.anneal()
print('最短总距离为：', e)
print('最佳路径为：', state)
# 使用 matplotlib 绘图
order_by_best = list(map(int, state)) + [int(state[0])]
plt.plot([cities[i][0] for i in order_by_best], [cities[i][1] for i in order_by_best], 'xb-')
for city, pos in cities.items():
plt.text(pos[0], pos[1], str(city))
plt.show()
点击查看代码
%% I. 清空环境变量
clear all
clc
%% II. 一元函数优化
x = 1:0.01:2;
y = sin(10*pi*x) ./ x;
figure
plot(x,y,'linewidth',1.5)
ylim([-1.5, 1.5])
xlabel('x')
ylabel('y')
title('y = sin(10*pi*x) / x')
hold on
%%
% 1. 标记出最大值点
[maxVal,maxIndex] = max(y);
plot(x(maxIndex), maxVal, 'r*','linewidth',2)
text(x(maxIndex), maxVal, {['X: ' num2str(x(maxIndex))];['Y: ' num2str(maxVal)]})
hold on
%%
% 2. 标记出最小值点
[minVal,minIndex] = min(y);
plot(x(minIndex), minVal, 'ks','linewidth',2)
text(x(minIndex), minVal, {['X: ' num2str(x(minIndex))];['Y: ' num2str(minVal)]})
function fitnessVal = fitness( x )
fitnessVal = sin(10*pi*x) / x;
end
fun = @fitness;
x0 = [1];
lb = [1];
ub = [2];
x = simulannealbnd(fun,x0,lb,ub)
参考网址
https://www.zhihu.com/search?type=content&q=数学建模模拟退火
https://zhuanlan.zhihu.com/p/617855897
4 聚类分析
可能需要用到spss
\(K\)均值聚类算法
将样本划分为由类似对象组成的多个类，\(K\)值为类的个数，选取\(K\)个数据对象作为初始聚类中心，计算其余各个数据到各个中心的距离，并归到距离最短的那个类中
更新聚类中心（不一定是具体的数据点，可以是空的地方，类似于物理上的质心）
在初始化\(K\)个聚类中心时，可以用\(k means++\)进行优化，选取相距尽可能远的聚类中心
系统（层次）聚类
将每个对象看作一个类，所以初始一共有 \(n\) 个类
将距离最近的两个类合成一个新的类
不断重复，直到最后合为一个类
通过肘部法则判断最适合的\(K\)值
DBSCAN算法
参考网址
https://zhuanlan.zhihu.com/p/673950122
5 差分进化模型
6 决策树
7 灰色预测模型
解决问题
处理少量数据或缺乏历史数据，短期预测和趋势分析
基本原理
1.通过累加将原始数据序列转变为新数据序列
2.利用微分方程、矩阵求解预测值
3.利用方差进行精度检验
参考代码
点击查看代码
function []=huidu()
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法
y=input("请输入数据);
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
yy(i)=yy(i-1)+y(i);
end
for i=1:(n-1)
B(i,1)=-(yy(i)+yy(i+1))/2;B(i,2)=1;
B=ones(n-1,2);
end
BT=B';
for j=1:n-1
YN(j)=y(j+1);
end
YN=YN';
A=inv(BT*B)*BT*YN;
a=A(1);
U=A(2);
t=u/a;
i=l:n+2;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;yys(1)=y(1);
for j=n+2:-1:2
ys(j)=yys(j)-yys(j-1);
end
X=1:n;
xs=2:n+2;yn=ys(2:n+2);plot(x,y,'^r',xs,yn,'*-b');det=0;
sum1=0;
sumpe=0;
for i=l:n
sumpe=sumpe+y(i);
end
pe=sumpe/n;
for i=1:n
sum1=sum1+(y(i)-pe).^2;
end
s1=sqrt(sum1/n);
sumce=0;
for i=2:n
sumce=sumce+(y(i)-yn(i));
end
ce=sumce/(n-1);sum2=0;
for i=2:n
sum2=sum2+(y(i)-yn(i)-ce).^2;
end
s2=sqrt(sum2/(n-1));
c=(s2)/(s1);
disp(['后验差比值为:,num2str(c)]);
if c<8.35
disp(系统预测精度好 )
else if c<0.5
disp("系统预测精度合格 )
else if c<0.65
disp("系统预测精度勉强')
else
disp("系统预测精度不合格 )
end
end
end
disp(['下个拟合值为,num2str(ys(n+1))]);
disp([再下个拟合值为',num2str(ys(n+2))])
参考网址
https://zhuanlan.zhihu.com/p/637849579
浙公网安备 33010602011771号