数学分析模型总结

层次分析法

1.建立层次结构模型

将决策的目标、考虑的因素（决策准则）和决策对象按照他们之间的相互关系分为最高层、中间层和最低层，绘出层次结构图。 最高层： 决策的目的、要解决的问题。 最低层： 决策时的备选方案。 中间层： 考虑的因素、决策的准则。 对相邻的两层，称高层为目标层，低层为因素层。 层次分析法所要解决的问题是关于最低层对最高层的相对权重的问题，按此相对权重可以对最低层中的各种方案、措施进行排序，从而在不同的方案中做出选择或形成选择方案的原则。

2.构造判断矩阵 层次分析法中构造判断矩阵的方法是一致矩阵法，即：不把所有因素放在一起比较，而是两两相互比较；对此时采用相对尺度，以尽可能减少性质不同因素相互比较的困难，以提高准确度。

3.层次单排序及其一致性检验 对应于判断矩阵最大特征根λ m a x的特征向量，经归一化（使向量中各元素之和为1）后记为W。W的元素为同一层次元素对于上一层因素某因素相对重要性的排序权值，这一过程称为层次单排序。

4.层次总排序及其一致性检验

  • 计算某一层次所有因素对于最高层（总目标）相对重要性的权值，称为层次总排序。
  • 这一过程是从最高层次到最低层次依次进行的。
image-20240130105508933
image-20240130105522533
image-20240130105530952
image-20240130105602627
image-20240130105612247
image-20240130105622706
image-20240130105634826
image-20240130105642063

多属性决策模型

加权平均

image-20240130110100685

正向化标准化（内容不贴了）

image-20240130110331660
image-20240130110339299
image-20240130110351212
image-20240130110356088

灰色预测

美赛不推荐使用，评委不知道这个算法

最短路径算法

Dijkstra

image-20240130110629888

每次从 「未求出最短路径的点」中 取出 距离距离起点 最小路径的点，以这个点为桥梁 刷新「未求出最短路径的点」的距image-20240130110726035

代码实例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include<bits/stdc++.h>
using namespace std;
#define PII pair<int ,int>
const int INF=1e10;
const int N=1000010;
int n,m,s;
int dis[N];
int vis[N];
vector<PII>E[N];
void dj(int s)
{
	for(int i=1;i<=n;i++)
	{
		dis[i]=INF;
		vis[i]=0;
	}
	dis[s]=0;
	priority_queue<PII,vector<PII>,greater<PII> >q;
	q.push({0,s});
    while(!q.empty())
{
int t=q.top().second;
	q.pop();
		if(vis[t]){
			continue;
		}
		vis[t]=1;
	for(int i=0,l=E[t].size();i<l;i++)
	{
		int v=E[t][i].first;
		int w=E[t][i].second;
	
		if(dis[v]>dis[t]+w)
		{		
			dis[v]=dis[t]+w;
			q.push({dis[v],v});		
		}	
		}	
	}	
}

int main()
{
	 ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
cin>>n>>m>>s;
int u,v,w;
for(int i=0;i<m;i++)
{
	cin>>u>>v>>w;
	E[u].push_back({v,w});
	}	
	dj(s);
	for(int i=1;i<=n;i++)
	{
		cout<<dis[i]<<" ";
	}
	return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
%% Matlab作无向图
% （1）无权重（每条边的权重默认为1）
% 函数graph(s,t)：可在 s 和 t 中的对应节点之间创建边，并生成一个图
% s 和 t 都必须具有相同的元素数；这些节点必须都是从1开始的正整数，或都是字符串元胞数组。
s1 = [1,2,3,4];
t1 = [2,3,1,1];
G1 = graph(s1, t1);
plot(G1)
% 注意哦，编号最好是从1开始连续编号，不要自己随便定义编号
s1 = [1,2,3,4];
t1 = [2,3,1,1];
G1 = graph(s1, t1);
plot(G1)

% 注意字符串元胞数组是用大括号包起来的哦
s2 = {'学校','电影院','网吧','酒店'};
t2 = {'电影院','酒店','酒店','KTV'};
G2 = graph(s2, t2);
plot(G2, 'linewidth', 2)  % 设置线的宽度
% 下面的命令是在画图后不显示坐标
set( gca, 'XTick', [], 'YTick', [] );  

% （2）有权重
% 函数graph(s,t,w)：可在 s 和 t 中的对应节点之间以w的权重创建边，并生成一个图
s = [1,2,3,4];
t = [2,3,1,1];
w = [3,8,9,2];
G = graph(s, t, w);
plot(G, 'EdgeLabel', G.Edges.Weight, 'linewidth', 2) 
set( gca, 'XTick', [], 'YTick', [] );  

%% Matlab作有向图
% 无权图 digraph(s,t)
s = [1,2,3,4,1];
t = [2,3,1,1,4];
G = digraph(s, t);
plot(G)
set( gca, 'XTick', [], 'YTick', [] );  

% 有权图 digraph(s,t,w)
s = [1,2,3,4];
t = [2,3,1,1];
w = [3,8,9,2];
G = digraph(s, t, w);
plot(G, 'EdgeLabel', G.Edges.Weight, 'linewidth', 2) 
set( gca, 'XTick', [], 'YTick', [] );  

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
%% 注意：以下代码需要较新版本的matlab才能运行（最好是2016版本及以上哦）
% 如果运行出错请下载新版的matlab代码再运行

% 注意哦，Matlab中的图节点要从1开始编号，所以这里把0全部改为了9
% 编号最好是从1开始连续编号，不要自己随便定义编号
s = [9 9 1 1 2 2 2 7 7 6 6  5  5 4];
t = [1 7 7 2 8 3 5 8 6 8 5  3  4 3];
w = [4 8 3 8 2 7 4 1 6 6 2 14 10 9];
G = graph(s,t,w);
plot(G, 'EdgeLabel', G.Edges.Weight, 'linewidth', 2) 
set( gca, 'XTick', [], 'YTick', [] );  
[P,d] = shortestpath(G, 9, 4)  %注意：该函数matlab2015b之后才有哦

% 在图中高亮我们的最短路径
myplot = plot(G, 'EdgeLabel', G.Edges.Weight, 'linewidth', 2);  %首先将图赋给一个变量
highlight(myplot, P, 'EdgeColor', 'r')   %对这个变量即我们刚刚绘制的图形进行高亮处理（给边加上r红色）

% 求出任意两点的最短路径矩阵
D = distances(G)   %注意：该函数matlab2015b之后才有哦
D(1,2)  % 1 -> 2的最短路径
D(9,4)  % 9 -> 4的最短路径

% 找出给定范围内的所有点  nearest(G,s,d)
% 返回图形 G 中与节点 s 的距离在 d 之内的所有节点
[nodeIDs,dist] = nearest(G, 2, 10)   %注意：该函数matlab2016a之后才有哦

Floyd

image-20240130111841434

Floyd算法代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
function [dist,path] = Floyd_algorithm(D)
%% 该函数用于求解一个权重邻接矩阵任意两个节点之间的最短路径
% 输入：
%        D是权重邻接矩阵
% 输出：
%        dist是最短距离矩阵，其元素dist_ij表示表示i,j两个节点的最短距离
%        path是路径矩阵，其元素path_ij表示起点为i，终点为j的两个节点之间的最短路径要经过的节点

n = size(D,1);  % 计算节点的个数

% 初始化dist矩阵
dist = D;

% 下面我们来初始化path矩阵
path = zeros(n);
for j = 1:n
    path(:,j) = j;   % 将第j列的元素变为j
end
for i = 1:n
    path(i,i) = -1;  % 将主对角线元素变为-1
end

% 下面开始三个循环
for k=1:n    % 中间节点k从1- n 循环
   for i=1:n     % 起始节点i从1- n 循环
      for j=1:n    % 终点节点j从1-n 循环
          if dist(i,j)>dist(i,k)+dist(k,j)  % 如果i,j两个节点间的最短距离大于i和k的最短距离+k和j的最短距离
             dist(i,j)=dist(i,k)+dist(k,j);  % 那么我们就令这两个较短的距离之和取代i,j两点之间的最短距离
             path(i,j)=path(i,k);   % 起点为i，终点为j的两个节点之间的最短路径要经过的节点更新为path(i,k)
             % 注意，上面一行语句不能写成path(i,j) = k; 这是网上很多地方都容易犯的错误，在PPT11页中会告诉大家为什么不能这么写
          end
      end
   end
end

end

例子

有向图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
% PPT第七页的例子
%% 首先将图转换为权重邻接矩阵D
n = 5;  %一共五个节点
D = ones(n) ./ zeros(n);  % 全部元素初始化为Inf
for i = 1:n
    D(i,i) = 0;  % 主对角线元素为0
end
D(1,2) = 3;
D(1,3) = 8;
D(1,5) = -4;
D(2,5) = 7;
D(2,4) = 1;
D(3,2) = 4;
D(4,3) = -5;
D(5,4) = 6;
D(4,1) = 2;

%% 调用Floyd_algorithm函数求解
[dist,path] = Floyd_algorithm(D)

print_path(path,dist,1,5)
print_path(path,dist,1,4)
print_path(path,dist,3,1)

clc
disp('下面我们打印任意两点之间的最短距离：')
print_all_path(D)

无向图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
% 思考题的参考答案
%% 首先将图转换为权重邻接矩阵D
n = 9;  %一共九个节点
D = zeros(n);  % 全部元素初始化为0, 等会你们就知道为什么这样设置啦
% 因为是无向图，所以权重邻接矩阵是一个对称矩阵
D(1,2) = 4; D(1,8) = 8; 
D(2,8) = 3; D(2,3) = 8;
D(8,9) = 1; D(8,7) = 6; 
D(9,7) = 6; D(9,3) = 2;
D(7,6) = 2; D(3,4) = 7; 
D(3,6) = 4; D(6,4) = 14;
D(4,5) = 9; D(6,5) = 10;
D = D+D';   % 这个操作可以得到对称矩阵的另一半
for i = 1:n
    for j = 1:n
        if (i ~= j) && (D(i,j) == 0)  
            D(i,j) = Inf;   % 将非主对角线上的0元素全部变为Inf
        end
    end
end

%% 调用Floyd_algorithm函数求解
[dist,path] = Floyd_algorithm(D)


打印

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function [] = print_all_path(D)
%% 该函数的作用是求解一个权重邻接矩阵任意两个节点之间的最短路径,并打印所有的结果出来
% 输入：
%        D是权重邻接矩阵
% 输出：无

[dist,path] = Floyd_algorithm(D);   % 调用之前的Floyd_algorithm函数
n = size(D,1);
if n == 1
    warning('请输入至少两阶以上的权重邻接矩阵')   % 在屏幕中提示警告信息
    return;   % 不运行下面的语句，直接退出函数
end

for i = 1:n
    for j = 1:n
        if i ~= j  % 不等号用~=表示
            print_path(path,dist,i,j);   % 调用之前的print_path函数
            disp('-------------------------------------------')
            disp('  ')
        end
    end
end

end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
function [] = print_path(path,dist,i,j)
%% 该函数的作用是打印从i到j经过的最短路径
% 输入：
%        path是使用floyd算法求出来的路径矩阵
%        dist是使用floyd算法求出来的最短距离矩阵
%        i是起始节点的编号
%        j是终点节点的编号
% 输出：无

if i == j
    warning('起点和终点相同，请检查后重新输入')  % 在屏幕中提示警告信息
    return;  % 不运行下面的语句，直接退出函数
end
if path(i,j) == j   % 如果path(i,j) = j，则有两种可能：
% （1）如果dist(i,j) 为 Inf , 则说明从i到j没有路径可以到达
    if dist(i,j) == Inf
        disp(['从',num2str(i),'到',num2str(j),'没有路径可以到达'])
% （2）如果dist(i,j) 不为 Inf , 则说明从i到j可直接到达，且为最短路径
    else
        disp(['从',num2str(i),'到',num2str(j),'的最短路径为'])
        disp([num2str(i),' ---> ',num2str(j)])
        disp(['最短距离为',num2str(dist(i,j))])
    end
else  % 如果path(i,j) ~= j，则说明中间经过了其他节点：
    k = path(i,j);
    result = [num2str(i),' ---> '];  % 初始化要打印的这个字符串
    while k ~= j  % 只要k不等于j, 就一直循环下去
        result = [result , num2str(k) , ' ---> ' ];  % i先走到k这个节点处
        k = path(k,j);
    end
    result = [result , num2str(k)];
    disp(['从',num2str(i),'到',num2str(j),'的最短路径为'])
    disp(result)
    disp(['最短距离为',num2str(dist(i,j))])
end

end


模拟退火算法

模拟退火是物理上退火的方法，通过N次迭代，逼近函数上的一个值

大方向：循环算法

模拟退火算法(Simulated Annealing，简称SA)的思想最早是由Metropolis等提出的。其出发点是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性。模拟退火法是一种通用的优化算法，其物理退火过程由以下三部分组成:

（1） 加温过程。其目的是增强粒子的热运动，使其偏离平衡位置。当温度足够高时，固体将熔为液体，从而消除系统原先存在的非均匀状态。

（2） 等温过程。对于与周围环境交换热量而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行的，当自由能达到最小时，系统达到平衡状态。

（3） 冷却过程。使粒子热运动减弱，系统能量下降，得到晶体结构。

加温过程相当于对算法设定初值，等温过程对应算法的Metropolis抽样过程，冷却过程对应控制参数的下降。这里能量的变化就是目标函数，我们要得到的最优解就是能量最低态。其中Metropolis准则是SA算法收敛于全局最优解的关键所在，Metropolis准则以一定的概率接受恶化解，这样就使算法跳离局部最优的陷阱。

SA算法的Metropolis准则允许接受一定的恶化解，具体来讲，是以一定概率来接受非最优解。举个例子，相当于保留一些“潜力股”，使解空间里有更多的可能性。对比轮盘赌法，从概率论来讲，它是对非最优解给予概率0，即全部抛弃。

模拟退火本身是求一个最小值问题，但可以转化为求最大值问题，只需要对目标函数加个负号或者取倒数。

img点击并拖拽以移动编辑

不断滚动，概率变动

imgimg

img 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
%% SA 模拟退火: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值(动画演示）
tic
clear; clc
 
%% 绘制函数的图形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
figure
plot(x,y,'b-')
hold on  % 不关闭图形，继续在上面画图
 
%% 参数初始化
narvs = 1; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 200;  % 最大迭代次数
Lk = 100;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数
x_lb = -3; % x的下界
x_ub = 3; % x的上界
 
%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
y0 = Obj_fun1(x0); % 计算当前解的函数值
h = scatter(x0,y0,'*r');  % scatter是绘制二维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）
 
%% 定义一些保存中间过程的量，方便输出结果和画图
max_y = y0;     % 初始化找到的最佳的解对应的函数值为y0
MAXY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_y (方便画图）
 
%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  % 内循环，在每个温度下开始迭代
        y = randn(1,narvs);  % 生成1行narvs列的N(0,1)随机数
        z = y / sqrt(sum(y.^2)); % 根据新解的产生规则计算z
        x_new = x0 + z*T; % 根据新解的产生规则计算x_new的值
        % 如果这个新解的位置超出了定义域，就对其进行调整
        for j = 1: narvs
            if x_new(j) < x_lb(j)
                r = rand(1);
                x_new(j) = r*x_lb(j)+(1-r)*x0(j);
            elseif x_new(j) > x_ub(j)
                r = rand(1);
                x_new(j) = r*x_ub(j)+(1-r)*x0(j);
            end
        end
        x1 = x_new;    % 将调整后的x_new赋值给新解x1
        y1 = Obj_fun1(x1);  % 计算新解的函数值
        if y1 > y0    % 如果新解函数值大于当前解的函数值
            x0 = x1; % 更新当前解为新解
            y0 = y1;
        else
            p = exp(-(y0 - y1)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                x0 = x1; % 更新当前解为新解
                y0 = y1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if y0 > max_y  % 如果当前解更好，则对其进行更新
            max_y = y0;  % 更新最大的y
            best_x = x0;  % 更新找到的最好的x
        end
    end
    MAXY(iter) = max_y; % 保存本轮外循环结束后找到的最大的y
    T = alfa*T;   % 温度下降
    pause(0.01)  % 暂停一段时间(单位：秒)后再接着画图
    h.XData = x0;  % 更新散点图句柄的x轴的数据（此时解的位置在图上发生了变化）
    h.YData = Obj_fun1(x0); % 更新散点图句柄的y轴的数据（此时解的位置在图上发生了变化）
end
 
disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(max_y)
 
pause(0.5)
h.XData = [];  h.YData = [];  % 将原来的散点删除
scatter(best_x,max_y,'*r');  % 在最大值处重新标上散点
title(['模拟退火找到的最大值为', num2str(max_y)])  % 加上图的标题
 
%% 画出每次迭代后找到的最大y的图形
figure
plot(1:maxgen,MAXY,'b-');
xlabel('迭代次数');
ylabel('y的值');
toc

旅行商问题也可以用蒙特卡罗解决

种群竞争模型

image-20240130112916180
image-20240130112922246 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
clc;clear
%  Matlab求不出来解析解
% dsolve('Dx1 = 0.5*x1*(1-x1/300-0.5*x2/500)','Dx2=0.5*x2*(1-x2/500-2*x1/300)','x1(0)=80,x2(0)=100','t')

% 下面用ode45函数求数值解
% 自变量为时间t，范围为0-30； 甲乙两个种群的数量初始值为80，100（随便给的，大家可以调整来看结果的变化）
[t,x]=ode45('fun',[0 30],[80 100]); 
plot(t,x(:,1),'r-',t,x(:,2),'b-')  % x的第一列是甲种群数量，x的第二列是乙种群数量
legend('种群甲','种群乙')
% axis([0 30 0 500])


function dx=fun(t,x)   % 大家可以修改里面的参数，来看结果的变化
    r1=0.5;  r 2=0.5; % 甲乙的增长率
%     r1=0.8;  r2=1; % 甲乙的增长率
    N1=300;   N2=500;   % 甲乙的最大数量
    % sigma1: 单位数量的乙种群(相对于N2)消耗的供养甲的食物量为单位数量的甲(相对于N1)消耗的供养甲的食物量的倍数。
    % sigma2: 单位数量的甲种群(相对于N1)消耗的供养乙的食物量为单位数量的乙(相对于N2)消耗的供养乙的食物量的倍数。
    sigma1=0.5;  sigma2=2;   
%     sigma1=0.5;  sigma2=4;   
%     sigma1=0.4;  sigma2=0.2;
%   当sigma1和sigma2同时大于1时（这种现象本身在自然界就几乎不可能出现），得到的结果不稳定。
%     sigma1=3;  sigma2=2;   
%     sigma1=2.2;  sigma2=2;   

    dx = zeros(2,1);
    dx(1) = r1*x(1)*(1-x(1)/N1-sigma1*x(2)/N2);
    dx(2) = r2*x(2)*(1-x(2)/N2-sigma2*x(1)/N1);
end

image-20240130113132533

数学规划模型

概述

image-20240128212145083
image-20240128212202777
image-20240128212217297
image-20240128212234617

线性规划

image-20240128212439349
image-20240128212611139
image-20240128212925755

代码命令解释

image-20240128213014243

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
%% Matlab求解线性规划
% [x fval] = linprog(c, A, b, Aeq, beq, lb,ub, x0)  
% c是目标函数的系数向量，A是不等式约束Ax<=b的系数矩阵，b是不等式约束Ax<=b的常数项
% Aeq是等式约束Aeq x=beq的系数矩阵，beq是等式约束Aeq x=beq的常数项
% lb是X的下限，ub是X的上限，X是向量[x1,x2,...xn]' , 即决策变量。
% 迭代的初始值为x0（一般不用给）
% 更多该函数的用法说明请看讲义

%% 例题1
c = [-5 -4 -6]';  % 加单引号表示转置
% c = [-5 -4 -6];  % 写成行向量也是可以的，不过不推荐，我们按照标准型来写看起来比较正规
A = [1 -1 1;
        3 2 4;
        3 2 0];
b = [20 42 30]';   
lb = [0 0 0]'; 
[x fval] = linprog(c, A, b, [], [], lb)  % ub我们直接不写，则意味着没有上界的约束
% x =
%          0
%    15.0000
%     3.0000
% 
% fval =
%    -78


%% 例题2
c = [0.04 0.15 0.1 0.125]';  
A = [-0.03 -0.3 0 -0.15;
        0.14 0 0 0.07];
b = [-32 42]';
Aeq = [0.05 0 0.2 0.1];
beq = 24;
lb = [0 0 0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)
% x =
%          0
%   106.6667
%   120.0000
%          0
% 
% fval =
%     28

% 这个题可能有多个解，即有多个x可以使得目标函数的最小值为28（不同的Matlab版本可能得到的x的值不同，但最后的最小值一定是28）
% 例如我们更改一个限定条件：令x1要大于0（注意Matlab中线性规划的标准型要求的不等式约束的符号是小于等于0）
% x1 >0  等价于  -x1 < 0，那么给定 -x1 <= -0.1 (根据实际问题可以给一个略小于0的数-0.1)，这样能将小于号转换为小于等于号，满足Matlab的标准型
c = [0.04 0.15 0.1 0.125]';  
A = [-0.03 -0.3 0 -0.15;
        0.14 0 0 0.07
        -1 0 0 0];
b = [-32 42 -0.1]';
Aeq = [0.05 0 0.2 0.1];
beq = 24;
lb = [0 0 0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)
% x =
%     0.1000
%   106.6567
%   119.9750
%          0
%
% fval =
%    28.0000


%% 例题3
c = [-2 -3 5]';
A = [-2 5 -1;
          1 3 1];
b = [-10 12];
Aeq = ones(1,3);
beq = 7;
lb = zeros(3,1);
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval % 注意这个fval要取负号（原来是求最大值，我们添加负号变成了最小值问题）
% x =
%     6.4286
%     0.5714
%          0
% fval =
%   -14.5714
% fval =
%    14.5714


%% 多个解的情况
% 例如 ： min z = x1 + x2   s.t.  x1 + x2 >= 10
c = [1 1]';   
A = [-1 -1];
b = -10;
[x fval] = linprog(c, A, b)   % Aeq, beq, lb和ub我们都没写，意味着没有等式约束和上下界约束
% x有多个解时，Matlab会给我们返回其中的一个解

%% 不存在解的情况
% 例如 ： min z = x1 + x2   s.t.  x1 + x2 = 10 、 x1 + 2*x2 <= 8、 x1 >=0 ，x2 >=0 
c = [1 1]'; 
A = [1 2];
b = 8;
Aeq = [1 1];
beq = 10;
lb = [0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)  % Linprog stopped because no point satisfies the constraints.（没有任何一个点满足约束条件）

线性规划的典型例题

image-20240128213640923 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
%% 生产决策问题
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
% (1) 系数向量
c = zeros(9,1); % 初始化目标函数的系数向量全为0
c(1) = 1.25 -0.25 -300/6000*5;  % x1前面的系数是c1
c(2) = 1.25 -0.25 -321/10000*7;
c(3) = -250 / 4000 * 6;
c(4)  = -783/7000*4;
c(5) = -200/4000 * 7;
c(6) = -300/6000*10;
c(7) = -321 / 10000 * 9;
c(8) = 2-0.35-250/4000*8;
c(9) = 2.8-0.5-321/10000*12-783/7000*11;
c = -c;  % 我们求的是最大值，所以这里需要改变符号
% (2) 不等式约束
A = zeros(5,9);
A(1,1) = 5;  A(1,6) = 10;
A(2,2) = 7;  A(2,7) = 9; A(2,9) = 12;
A(3,3) = 6;  A(3,8) = 8;
A(4,4) = 4;  A(4,9) = 11;
A(5,5) = 7;  
b = [6000 10000 4000 7000 4000]';
% (3) 等式约束
Aeq = [1 1 -1 -1 -1 0 0 0 0;
            0 0 0 0 0 1 1 -1 0];
beq = [0 0]';
%（4）上下界
lb = zeros(9,1);

% 进行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval
% fval =
%           1146.56650246305
%  注意，本题应该是一个整数规划的例子，我们在后面的整数规划部分再来重新求解。
intcon = 1:9;
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb)
fval = -fval

整数规划

image-20240128214116824
image-20240128214417024 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
%% 线性整数规划问题
%% 例1
c=[-20,-10]';
intcon=[1,2];  % x1和x2限定为整数
A=[5,4;
      2,5];
b=[24;13];
lb=zeros(2,1);  
[x,fval]=intlinprog(c,intcon,A,b,[],[],lb)
fval = -fval

%% 例2
c=[18,23,5]';
intcon=3;  % x3限定为整数
A=[107,500,0;
      72,121,65;
      -107,-500,0;
      -72,-121,-65];
b=[50000;2250;-500;-2000];
lb=zeros(3,1);
[x,fval]=intlinprog(c,intcon,A,b,[],[],lb)

%% 例3
c=[-3;-2;-1]; intcon=3; % x3限定为整数
A=ones(1,3); b=7;
Aeq=[4 2 1]; beq=12;
lb=zeros(3,1); ub=[+inf;+inf;1]; %x(3)为0-1变量
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)


背包问题

image-20240128214732168 
1
2
3
4
5
6
7
8
9
10
11
12
13
%% 背包问题（货车运送货物的问题）
c = -[540 200 180 350 60 150 280 450 320 120];  % 目标函数的系数矩阵(最大化问题记得加负号)
intcon=[1:10];  % 整数变量的位置(一共10个决策变量，均为0-1整数变量)
A = [6 3 4 5 1 2 3 5 4 2];  b = 30;   % 线性不等式约束的系数矩阵和常数项向量（物品的重量不能超过30）
Aeq = []; beq =[];  % 不存在线性等式约束
lb = zeros(10,1);  % 约束变量的范围下限
ub = ones(10,1);  % 约束变量的范围上限
%最后调用intlinprog()函数
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
fval = -fval



指派问题

image-20240128215024353 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
%% 指派问题（选择队员去进行游泳接力比赛）
clear;clc
c = [66.8 75.6 87 58.6 57.2 66 66.4 53 78 67.8 84.6 59.4 70 74.2 69.6 57.2 67.4 71 83.8 62.4]';  % 目标函数的系数矩阵（先列后行的写法）
intcon = [1:20];  % 整数变量的位置(一共20个决策变量，均为0-1整数变量)
% 线性不等式约束的系数矩阵和常数项向量（每个人只能入选四种泳姿之一，一共五个约束）
A = [1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1];
% A = zeros(5,20);
% for i = 1:5
%     A(i, (4*i-3): 4*i) = 1;
% end
b = [1;1;1;1;1];
% 线性等式约束的系数矩阵和常数项向量 （每种泳姿有且仅有一人参加，一共四个约束）
Aeq = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;
          0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0;
          0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0;
          0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1];
% Aeq = [eye(4),eye(4),eye(4),eye(4),eye(4)];  % 或者写成 repmat(eye(4),1,5)  
beq = [1;1;1;1];
lb = zeros(20,1);  % 约束变量的范围下限
ub = ones(20,1);  % 约束变量的范围上限
%最后调用intlinprog()函数
[x,fval] = intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
% reshape(x,4,5)'
%      0     0     0     1    甲自由泳
%      1     0     0     0    乙蝶泳
%      0     1     0     0    丙仰泳
%      0     0     1     0    丁蛙泳
%      0     0     0     0    戊不参加


钢管切割问题

image-20240128215414471
image-20240128215510504 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%% 钢管切割问题
%% (1)枚举法找出同一个原材料上所有的切割方法
for i = 0: 2  % 2.9m长的圆钢的数量
    for j = 0: 3  % 2.1m长的圆钢的数量
        for k = 0:6   % 1m长的圆钢的数量
            if 2.9*i+2.1*j+1*k >= 6 && 2.9*i+2.1*j+1*k <= 6.9
                disp([i, j, k])
            end
        end
    end
end
% 有同学使用比较老的MATLAB版本，会出现浮点数计算的误差
% 只需要将上面的if这一行进行适当的放缩即可。
% if 2.9*i+2.1*j+1*k >= 6-0.0000001 && 2.9*i+2.1*j+1*k <= 6.9+0.0000001
% 有兴趣的同学可以百度下：浮点数计算误差

%% (2) 线性整数规划问题的求解
c = ones(7,1);  % 目标函数的系数矩阵
intcon=[1:7];  %  整数变量的位置(一共7个决策变量，均为整数变量)
A = -[1 2 0 0 0 0 1;  
         0 0 3 2 1 0 1;
         4 1 0 2 4 6 1];  % 线性不等式约束的系数矩阵
b = -[100 100 100]'; %  线性不等式约束的常数项向量
lb = zeros(7,1); % 约束变量的范围下限
[x,fval]=intlinprog(c,intcon,A,b,[],[],lb)



非线性问题的求解

image-20240128215640424
image-20240128215838554
image-20240128220006524

代码求解

可以先给定不同初始值，在里面找到最优解

也可以蒙特卡罗模拟找到一个蒙特卡罗解，再作为初始值进行求解。

min最小值

f函数

con约束

X0是一个初始值，线性规划里边初始值对于结果没影响，而非线性规划中x0的选取很关键，因为求出的是一个局部最优解。

image-20240128221553606

求解的方法有四种，可以提高结果的稳健性能。

image-20240128221841171

例题一

如果考察了，建议蒙特卡洛加上四种方法一起用，这样就稳啦

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
%% 非线性规划的函数
% [x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
% x0表示给定的初始值（用行向量或者列向量表示），必须得写
% A b表示线性不等式约束
% Aeq beq 表示线性等式约束
% lb ub 表示上下界约束
% @fun表示目标函数
% @nonlfun表示非线性约束的函数
% option 表示求解非线性规划使用的方法
clear;clc
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）

%% 例题1的求解
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
% s.t. -(x1-1)^2 +x2 >= 0 ;  2x1-3x2+6 >= 0
x0 = [0 0];  %任意给定一个初始值 
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)  % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
fval = -fval
% 一个值得讨论的地方，能不能把线性不等式约束Ax <= b也写到nonlfun1函数中？
% 先把nonlfun1中的c改为下面这样：
% c = [(x(1)-1)^2-x(2); 
%        -2*x(1)+3*x(2)-6];
%  [x,fval] = fmincon(@fun1,x0,[],[],[],[],[],[],@nonlfun1)
% 结果也是可以计算出来的，但并不推荐这样做~

目标函数
function f = fun1(x)
    % 注意：这里的f实际上就是目标函数，函数的返回值也是f
    % 输入值x实际上就是决策变量，由x1和x2组成的向量
    % fun1是函数名称，到时候会被fmincon函数调用, 可以任意取名
    % 保存的m文件和函数名称得一致，也要为fun1.m
%      max  f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
    f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ; 
end
这是非线性约束

function [c,ceq] = nonlfun1(x)
    % 注意：这里的c实际上就是非线性不等式约束，ceq实际上就是非线性等式约束
    % 输入值x实际上就是决策变量，由x1和x2组成的一个向量
    % 返回值有两个，一个是非线性不等式约束c，一个是非线性等式约束ceq
    % nonlfun1是函数名称，到时候会被fmincon函数调用, 可以任意取名，但不能和目标函数fun1重名
    % 保存的m文件和函数名称得一致，也要为nonlfun1.m
%     -(x1-1)^2 +x2 >= 0 
   c = [(x(1)-1)^2-x(2)];   % 千万別写成了: (x1-1)^2 -x2
   ceq = [];  % 不存在非线性等式约束，所以用[]表示
end




%% 使用其他算法对例题1求解
% edit fmincon  % 查看fmincon的“源代码”
% Matlab2017a默认使用的算法是'interior-point' 内点法
% 使用interior point算法 （内点法）
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval
% 使用SQP算法 （序列二次规划法）
option = optimoptions('fmincon','Algorithm','sqp')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval   %得到-4.358，远远大于内点法得到的-1,猜想是初始值的影响
% 改变初始值试试
x0 = [1 1];  %任意给定一个初始值 
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  % 最小值为-1，和内点法相同（这说明内点法的适应性要好）
fval = -fval  
% 使用active set算法 （有效集法）
option = optimoptions('fmincon','Algorithm','active-set')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval  
% 使用trust region reflective (信赖域反射算法)
option = optimoptions('fmincon','Algorithm','trust-region-reflective')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval  
% this algorithm does not solve problems with the constraints you have specified. 
% 这说明这个算法不适用我们这个约束条件，所以以后遇到了不能求解的情况，记得更换其他算法试试！！！

%% 选取初始值得到的结果可能会不满足限定条件，出现了一个Bug 因此选择的初始值很重要
x0 = [40.8, 10.8];
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval  
% https://cn.mathworks.com/help/optim/ug/fmincon.html

%% 生成不同的随机初始值来优化代码，有一定几率会触发上面那个Bug，因此不推荐
n = 10;  % 重复n次
Fval = +inf; X = [0,0];  %初始化最优的结果
A = [-2 3]; b = 6;
for i = 1:n
    x0 = [rand()*10 , rand()*10];  %用随机数生成一个初始值（随机数的范围自己根据题目条件设置） 
    [x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option); % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
    if fval < Fval  % 如果找到了更小的值，那么就代替最优的结果
        Fval = fval;
        X = x;
    end
end
Fval = -Fval
X

%% 使用蒙特卡罗的方法来找初始值(推荐）
clc,clear;
n=10000000; %生成的随机数组数
x1=unifrnd(-100,100,n,1);  % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(-100,100,n,1);  % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x2
fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
for i=1:n
    x = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2]
    if ((x(1)-1)^2-x(2)<=0)  & (-2*x(1)+3*x(2)-6 <= 0)     % 判断是否满足条件
        result = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;  % 如果满足条件就计算函数值
        if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
            fmin = result;  % 那么就更新这个函数值为新的最小值
            x0 = x;  % 并且将此时的x1 x2更新为初始值
        end
    end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
fval = -fval

例题二

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
%% 例题二的求解
x0 = [1 1 1];  %任意给定一个初始值 
lb = [0 0 0];  % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
% x =
%          0.552167405729277          1.20325915507969         0.947824046150443
% fval =
%           10.6510918606939



%% 使用蒙特卡罗的方法来找初始值(推荐）
clc,clear;
n=1000000; %生成的随机数组数
x1= unifrnd(0,2,n,1);   % 生成在[0,2]之间均匀分布的随机数组成的n行1列的向量构成x1
x2 = sqrt(2-x1);  % 根据非线性等式约束用x1计算出x2
x3 = sqrt((3-x2)/2); % 根据非线性等式约束用x2计算出x3
fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
for i=1:n
    x = [x1(i), x2(i), x3(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
    if (-x(1)^2+x(2)-x(3)^2<=0) & (x(1)+x(2)^2+x(3)^2-20<=0)   % 判断是否满足条件
        result =sum(x.*x) + 8 ;  % 如果满足条件就计算函数值
        if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
            fmin = result;  % 那么就更新这个函数值为新的最小值
            x0 = x;  % 并且将此时的x1 x2 x3更新为初始值
        end
    end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)
lb = [0 0 0];  % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下





function f = fun2(x)
    %     f = x(1)^2+x(2)^2 +x(3)^2+8 ; 
    f = sum(x.*x) + 8;  % 可别忘了x实际上是一个向量，我们可以使用矩阵的运算符号对其计算
end

例题三

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
%% 例题三的求解(蒙特卡罗模拟那一讲的例题)
clear;clc
% 蒙特卡罗模拟得到的最大值为3445.6014
% 最大值处x1 x2 x3的取值为：
%           22.5823101903968          12.5823101903968          12.1265223966757
A = [1 -2 -2;  1 2 2];  b = [0 72];
x0 = [ 22.58   12.58  12.13];
Aeq = [1 -1 0]; beq = 10;
lb = [-inf 10 -inf];  ub = [inf 20 inf];  
[x,fval] = fmincon(@fun3,x0,A,b,Aeq,beq,lb,ub,[])  % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
fval = -fval


function f = fun3(x)
    f = -prod(x);  % 可别忘了x实际上是一个向量（prod表示连乘符号，用法和sum类似）
end



选址问题

背景

image-20240129110524287
image-20240129110713468 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
%% 选址问题
clear;clc
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
% % (1) 系数向量（原来线性规划问题的写法,我们只需要在此基础上改动一点就可以了）
% a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
% b=[1.25  0.75  4.75	5  6.5  7.25];   % 工地的纵坐标
% x = [5  2];  % 料场的横坐标
% y = [1  7];  % 料场的纵坐标
% c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
% for  j =1:2
%     for i = 1:6
%         c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
%     end
% end
% (2) 不等式约束
A =zeros(2,16);  % 注意这里要改成16
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式约束
Aeq = zeros(6,16);  % 注意这里要改成16
for i = 1:6
    Aeq(i,i) = 1;  Aeq(i,i+6) = 1;
end
beq = [3 5 4 7 6 11]';  % 每个工地的日需求量
%（4）上下界
lb = zeros(16,1);
% lb = [zeros(12,1); -inf*ones(4,1)];  两个新料场坐标的下界可以设为-inf

% 进行求解
% 注意哦，这里我们只尝试了这一个初始值，大家可以试试其他的初始值，有可能能够找到更好的解。
% 未来我会在遗传算法中再来看这个例题。
x0 = [3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7];  % 用第一问的结果作为初始值
[x,fval] = fmincon(@fun5,x0,A,b,Aeq,beq,lb)  % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
reshape(x(1:12),6,2)  % 将x的前12个元素变为6行2列便于观察（reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1）
% 新坐标（5.74，4.99） （7.25，7.25）
% fval =
%           89.9231692432933
% 第一问的fval =
%           135.281541790676
135.281541790676 - 89.9231692432933  %  45.3583725473827



function f = fun5(xx)  % 注意为了避免和下面的x同号，我们把决策变量的向量符号用xx表示（注意xx的长度为16）
    a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
    b=[1.25  0.75  4.75	5  6.5  7.25];   % 工地的纵坐标
    x = [xx(13)  xx(15)];  % 新料场的横坐标
    y = [xx(14)  xx(16)];  % 新料场的纵坐标
    c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
    for  j =1:2
        for i = 1:6
            c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
        end
    end
    % 下面我们要求吨千米数，注意c是列向量，我们计算非线性规划时给定的初始值x0是行向量
    f = xx(1:12) * c;
end







最大最小化模型

image-20240129005314765

典型例题

image-20240129005430279

模型求解

image-20240129005556115 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
%% 最大最小化模型  :   min{max[f1,f2,···,fm]}
x0 = [6, 6];      % 给定初始值
lb = [3, 4];  % 决策变量的下界
ub = [8, 10];  % 决策变量的上界
[x,feval] = fminimax(@Fun,x0,[],[],[],[],lb,ub)
max(feval)
% x =
%     8.0000    8.5000
% feval =
%    13.5000    5.5000    5.5000   12.5000    8.5000    8.5000    5.5000   13.5000    9.5000    0.5000
% 结论：
% 在坐标为(8,8.5)处建立供应中心可以使该点到各需求点的最大距离最小，最小的最大距离为13.5单位。



function f = Fun(x)
    a=[1 4 3 5 9 12 6 20 17 8];
    b=[2 10 8 18 1 4 5 10 8 9];
    %  函数向量
    f=zeros(10,1);
    for i = 1:10
        f(i) = abs(x(1)-a(i))+abs(x(2)-b(i));  
    end
% f(1) = abs(x(1)-a(1))+abs(x(2)-b(1));  
% f(2) = abs(x(1)-a(2))+abs(x(2)-b(2));
% f(3) = abs(x(1)-a(3))+abs(x(2)-b(3));
% f(4) = abs(x(1)-a(4))+abs(x(2)-b(4));
% f(5) = abs(x(1)-a(5))+abs(x(2)-b(5));
% f(6) = abs(x(1)-a(6))+abs(x(2)-b(6));
% f(7) = abs(x(1)-a(7))+abs(x(2)-b(7));
% f(8) = abs(x(1)-a(8))+abs(x(2)-b(8));
% f(9) = abs(x(1)-a(9))+abs(x(2)-b(9));
% f(10) = abs(x(1)-a(10))+abs(x(2)-b(10));
end

多目标规划问题

背景

image-20240129010114001

可能标准化，正向化，还有权重

例题：

image-20240129010157820
image-20240129010334935
image-20240129010535787 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
%%  多目标规划问题
w1 = 0.4;  w2 = 0.6;  % 两个目标函数的权重  x1 = 5  x2 = 2
w1 = 0.5;  w2 = 0.5;  % 两个目标函数的权重  x1 = 5  x2 = 2
w1 = 0.3;  w2 = 0.7;  % 两个目标函数的权重  x1 = 1  x2 = 6
c = [w1/30*2+w2/2*0.4 ;w1/30*5+w2/2*0.3];  % 线性规划目标函数的系数
A = [-1 -1];  b = -7; % 不等式约束
lb = [0 0]'; ub = [5 6]'; % 上下界
[x,fval] = linprog(c,A,b,[],[],lb,ub)
f1 = 2*x(1)+5*x(2)
f2 = 0.4*x(1) + 0.3*x(2)


%% 敏感性分析 
clear;clc
W1 = 0.1:0.001:0.5;  W2 = 1- W1;  
n =length(W1);
F1 = zeros(n,1);  F2 = zeros(n,1);   X1 = zeros(n,1);  X2 = zeros(n,1);   FVAL = zeros(n,1);
A = [-1 -1];  b = -7; % 不等式约束
lb = [0 0]; ub = [5 6]; % 上下界
for i = 1:n
    w1 = W1(i);  w2 = W2(i);
    c = [w1/30*2+w2/2*0.4 ;w1/30*5+w2/2*0.3];  % 线性规划目标函数的系数
    [x,fval] = linprog(c,A,b,[],[],lb,ub);
    F1(i) = 2*x(1)+5*x(2);
    F2(i) = 0.4*x(1) + 0.3*x(2);
    X1(i) = x(1);
    X2(i) = x(2);
    FVAL(i) = fval;
end

% 「Matlab」“LaTex字符汇总”讲解：https://blog.csdn.net/Robot_Starscream/article/details/89386748
% 在图上可以加上数据游标，按住Alt加鼠标左键可以设置多个数据游标出来。
figure(1) 
plot(W1,F1,W1,F2)
xlabel('f_{1}的权重') 
ylabel('f_{1}和f_{2}的取值')
legend('f_{1}','f_{2}')

figure(2)
plot(W1,X1,W1,X2)
xlabel('f_{1}的权重') 
ylabel('x_{1}和x_{2}的取值')
legend('x_{1}','x_{2}')

figure(3)
plot(W1,FVAL)  % 看起来是两个直线组合起来的下半部分
xlabel('f_{1}的权重') 
ylabel('综合指标的值')

主成分分析

image-20240130113809842
image-20240130113819544

主成分分析法

综述：数据降维的方法

可以用一种线性变换的思想去理解，比如二维的一条直线，我们可以通过变换，使得这一条直线落在x或y轴上，达到降维的效果。

去中心化（把坐标原点放在数据中心）

找坐标系，找到数据方差最大的方向，就是第一主成分。（如果第一主成分不足以表达，就考虑吧选取第二个）

为了有效反映原来信息，第一主成分和第二主成分的协方差为0.以此类推可以获得p个主成分。这些主成分是互不相关，是依次递减的。

累计方差贡献率大于百分之80就可以了，或者特征根大于1就可以了。

？根据线性代数的知识，我们需要一则伸缩，二则旋转，伸缩不是问题

旋转的矩阵R又从何而来？

即是协方差矩阵的特征向量就是R。

协方差定义=

img
image-20240129173904717
image-20240129173917094
image-20240129173940409
image-20240129173956777

例题：

image-20240129174252433
image-20240129174259805
image-20240129174305765 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
clear;clc
 load data1.mat   % 主成分聚类
%  load data2.mat   % 主成分回归

% 注意，这里可以对数据先进行描述性统计
% 描述性统计的内容见第5讲.相关系数
[n,p] = size(x);  % n是样本个数，p是指标个数

%% 第一步：对数据x标准化为X
X=zscore(x);   % matlab内置的标准化函数（x-mean(x)）/std(x)

%% 第二步：计算样本协方差矩阵
R = cov(X);

%% 注意：以上两步可合并为下面一步：直接计算样本相关系数矩阵
R = corrcoef(x);
disp('样本相关系数矩阵为：')
disp(R)

%% 第三步：计算R的特征值和特征向量
% 注意：R是半正定矩阵，所以其特征值不为负数
% R同时是对称矩阵，Matlab计算对称矩阵时，会将特征值按照从小到大排列哦
% eig函数的详解见第一讲层次分析法的视频
[V,D] = eig(R);  % V 特征向量矩阵  D 特征值构成的对角矩阵


%% 第四步：计算主成分贡献率和累计贡献率
lambda = diag(D);  % diag函数用于得到一个矩阵的主对角线元素值(返回的是列向量)
lambda = lambda(end:-1:1);  % 因为lambda向量是从小大到排序的，我们将其调个头
contribution_rate = lambda / sum(lambda);  % 计算贡献率
cum_contribution_rate = cumsum(lambda)/ sum(lambda);   % 计算累计贡献率  cumsum是求累加值的函数
disp('特征值为：')
disp(lambda')  % 转置为行向量，方便展示
disp('贡献率为：')
disp(contribution_rate')
disp('累计贡献率为：')
disp(cum_contribution_rate')
disp('与特征值对应的特征向量矩阵为：')
% 注意：这里的特征向量要和特征值一一对应，之前特征值相当于颠倒过来了，因此特征向量的各列需要颠倒过来
%  rot90函数可以使一个矩阵逆时针旋转90度，然后再转置，就可以实现将矩阵的列颠倒的效果
V=rot90(V)';
disp(V)


%% 计算我们所需要的主成分的值
m =input('请输入需要保存的主成分的个数:  ');
F = zeros(n,m);  %初始化保存主成分的矩阵（每一列是一个主成分）
for i = 1:m
    ai = V(:,i)';   % 将第i个特征向量取出，并转置为行向量
    Ai = repmat(ai,n,1);   % 将这个行向量重复n次，构成一个n*p的矩阵
    F(:, i) = sum(Ai .* X, 2);  % 注意，对标准化的数据求了权重后要计算每一行的和
end

%% (1)主成分聚类 ： 将主成分指标所在的F矩阵复制到Excel表格，然后再用Spss进行聚类
% 在Excel第一行输入指标名称（F1,F2, ..., Fm）
% 双击Matlab工作区的F,进入变量编辑中，然后复制里面的数据到Excel表格
% 导出数据之后，我们后续的分析就可以在Spss中进行。

%%（2）主成分回归：将x使用主成分得到主成分指标，并将y标准化，接着导出到Excel，然后再使用Stata回归
% Y = zscore(y);  % 一定要将y进行标准化哦~
% 在Excel第一行输入指标名称（Y,F1, F2, ..., Fm）
% 分别双击Matlab工作区的Y和F,进入变量编辑中，然后复制里面的数据到Excel表格
% 导出数据之后，我们后续的分析就可以在Stata中进行。

聚类分析

image-20240129194932249

聚类不知道类别

K-means聚类算法

image-20240129194952521

流程图：

image-20240129195121059

优缺点分析

image-20240129195132844

K-means++算法

image-20240129195226437

操作

image-20240129195330285
image-20240129195405248

系统（层次）聚类

image-20240129195509250

距离计算

image-20240129195530809
image-20240129200236007
image-20240129200246282
image-20240129200254274
image-20240129200321590
image-20240129200335096

注意

image-20240129200430566
image-20240129200448939

会生成聚类谱系图，以此判断选择几类

image-20240129200741109
image-20240129200821668

详细操作见下博客

SPSS操作(四)：系统聚类分析_聚类分析spss操作-CSDN博客

DBSCAN 算法

image-20240129201052622

基本概念

image-20240129201312072
image-20240129201355010

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
clc;
clear;
close all;

%% Load Data

load mydata;


%% Run DBSCAN Clustering Algorithm

epsilon=0.5;
MinPts=10;
IDX=DBSCAN(X,epsilon,MinPts);




function [IDX, isnoise]=DBSCAN(X,epsilon,MinPts)

    C=0;
    
    n=size(X,1);
    IDX=zeros(n,1);  % 初始化全部为0，即全部为噪音点
    
    D=pdist2(X,X);
    
    visited=false(n,1);
    isnoise=false(n,1);
    
    for i=1:n
        if ~visited(i)
            visited(i)=true;
            
            Neighbors=RegionQuery(i);
            if numel(Neighbors)<MinPts
                % X(i,:) is NOISE
                isnoise(i)=true;
            else
                C=C+1;
                ExpandCluster(i,Neighbors,C);
            end
            
        end
    
    end
    
    function ExpandCluster(i,Neighbors,C)
        IDX(i)=C;
        
        k = 1;
        while true
            j = Neighbors(k);
            
            if ~visited(j)
                visited(j)=true;
                Neighbors2=RegionQuery(j);
                if numel(Neighbors2)>=MinPts
                    Neighbors=[Neighbors Neighbors2];   %#ok
                end
            end
            if IDX(j)==0
                IDX(j)=C;
            end
            
            k = k + 1;
            if k > numel(Neighbors)
                break;
            end
        end
    end
    
    function Neighbors=RegionQuery(i)
        Neighbors=find(D(i,:)<=epsilon);
    end

end

当然，让我逐行解释这段MATLAB代码：

  1. clc; clear; close all;
    • clc：清除命令窗口。
    • clear：清除工作区中的所有变量。
    • close all：关闭所有打开的图形窗口。
  2. load mydata;
    • 从名为 'mydata' 的文件中加载数据到工作区。这里的假设是 'mydata' 包含一个表示数据点的变量 X。
  3. epsilon=0.5; MinPts=10;
    • 定义DBSCAN算法的参数，epsilon 是邻域半径，MinPts 是邻域内最小数据点数。
  4. IDX=DBSCAN(X,epsilon,MinPts);
    • 调用DBSCAN函数，对数据 X 进行密度聚类，返回聚类结果 IDX。
  5. function [IDX, isnoise]=DBSCAN(X,epsilon,MinPts)
    • 定义DBSCAN算法的主函数，接受输入参数 X、epsilon 和 MinPts。
  6. C=0;
    • 初始化聚类簇数为0。
  7. n=size(X,1); IDX=zeros(n,1);
    • 获取数据点数量 n，初始化聚类标签 IDX 全部为0，表示所有点都是噪音点。
  8. D=pdist2(X,X);
    • 计算数据点之间的距离矩阵 D。
  9. visited=false(n,1); isnoise=false(n,1);
    • 初始化用于标记是否访问过的向量 visited 和标记是否为噪音点的向量 isnoise。
  10. for i=1:n
    • 开始对每个数据点进行迭代。
  11. if ~visited(i)
    • 如果当前点未被访问过，则执行以下操作。
  12. visited(i)=true; Neighbors=RegionQuery(i);
    • 将当前点标记为已访问，然后找到与当前点在邻域内的点集合 Neighbors。
  13. if numel(Neighbors)<MinPts
    • 如果邻域内点的数量小于 MinPts，则将当前点标记为噪音点。
  14. else
    • 否则，执行以下聚类操作。
  15. C=C+1; ExpandCluster(i,Neighbors,C);
    • 增加聚类簇数，并进行扩展聚类操作。
  16. function ExpandCluster(i,Neighbors,C)
    • 定义扩展聚类的子函数，给定当前点、邻域内点集合和当前簇数。
  17. IDX(i)=C;
    • 将当前点标记为属于当前簇。
  18. while true
    • 进入循环，不断扩展聚类。
  19. j = Neighbors(k);
    • 取出邻域内的第 k 个点。
  20. if ~visited(j)
    • 如果该点未被访问过，则执行以下操作。
  21. visited(j)=true; Neighbors2=RegionQuery(j);
    • 将该点标记为已访问，然后找到与该点在邻域内的点集合 Neighbors2。
  22. if numel(Neighbors2)>=MinPts
    • 如果新邻域内的点数量大于等于 MinPts，则将新邻域内的点添加到原邻域中。
  23. Neighbors=[Neighbors Neighbors2];
    • 将新邻域内的点添加到原邻域中。
  24. end
    • 结束新邻域内点的处理。
  25. if IDX(j)==0
    • 如果该点尚未被分配到任何簇，则将其分配到当前簇。
  26. IDX(j)=C;
    • 将该点标记为属于当前簇。
  27. k = k + 1; if k > numel(Neighbors) break; end
    • 处理邻域内的下一个点，直到邻域内的所有点都被处理完。
  28. function Neighbors=RegionQuery(i)
    • 定义邻域查询的子函数，给定当前点的索引 i，返回在邻域内的点的索引集合。
  29. Neighbors=find(D(i,:)<=epsilon);
    • 根据距离矩阵，找到与当前点距离在 epsilon 以内的点。
  30. end
    • 结束邻域查询子函数。
  31. end
    • 结束主函数。

这样，整个代码就实现了DBSCAN聚类算法。

多元回归分析

回归分析是数据分析中最基础也是最重要的分析工具，绝大多数的数据分析问题，都可以使用回归的思想来解决。回归分析的任务就是通过研究自变量X和因变量Y的相关关系，尝试去解释Y的形成机制，进而达到通过X去预测Y的目的

回归分析:研究X和Y相关性的分析（相关性≠因果性）

常见的回归分析有：线性回归、0-1回归、定序回归、计数回归和生存回归，其划分的依据是因变量y的类型。

image-20240129214522174
image-20240129214533470

回归分析的作用

image-20240129214553725

分类

image-20240129214629247

数据的分类： 横截面数据：在某一时点收集的不同对象的数据。

本章节主要是多元线性回归

image-20240129214810466

一元线性回归

image-20240129214828532

线性是灵动的

image-20240129214850518

注意引入变量时候要多加考虑

image-20240129215013434
image-20240129215038603

外生性的要求

image-20240129215112238

什么时候取对数

取对数的好处： （1）减弱数据的异方差性（2）如果变量本身不符合正态分布，取

了对数后 可能渐近服从正态分布（ 3 ）模型形式的需要，让模型具有经济学意义。

image-20240129215205907
image-20240129215212290

虚拟变量的解释:

image-20240129215246712
image-20240129215300615

多分类的虚拟变量的设置：

image-20240129215355156
image-20240129215401967

/为了避免完全多重共线性的影响，引入虚拟变量的个数一般是分类数减//1,另外一个为对照组//。/**

含有交互项的自变量：

image-20240129215424047

回归实例：

image-20240129215502987

操作步骤：

stata软件：

第一步：导入数据

image-20240129215527268

第二步：数据描述性统计

image-20240129215541264
image-20240129215607055
image-20240129215619067
image-20240129215627382
image-20240129215637199

拟合优度R²较低怎么办：

image-20240129215643497

标准回归化系数：

image-20240129215700385
image-20240129215707976

使用OLS时，扰动项μ需要满足的条件：

image-20240129215902585
image-20240129215928314
image-20240129215935706
image-20240129215945766
image-20240129220020904
image-20240129220048086
image-20240129220104904
image-20240129220124650
image-20240129220136200

多重共线性：

image-20240129220153354
image-20240129220202916

处理方法：

image-20240129220239152
image-20240129220248065
image-20240129220300213
image-20240129220322288

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
// 按键盘上的PageUp可以使用上一次输入的代码（Matlab中是上箭头）
// 清除所有变量
clear
// 清屏 和 matlab的clc类似
cls 
// 导入数据（其实是我们直接在界面上粘贴过来的，我们用鼠标点界面导入更方便 本条请删除后再复制到论文中，如果评委老师看到了就知道这不是你写的了）
// import excel "C:/Users/hc_lzp/Desktop/数学建模视频录制/第7讲.多元回归分析/代码和例题数据/课堂中讲解的奶粉数据.xlsx", sheet("Sheet1") firstrow
import excel "课堂中讲解的奶粉数据.xlsx", sheet("Sheet1") firstrow
// 定量变量的描述性统计
summarize 团购价元 评价量 商品毛重kg
// 定性变量的频数分布，并得到相应字母开头的虚拟变量
tabulate 配方,gen(A)
tabulate 奶源产地 ,gen(B)
tabulate 国产或进口 ,gen(C)
tabulate 适用年龄岁 ,gen(D)
tabulate 包装单位 ,gen(E)
tabulate 分类 ,gen(F)
tabulate 段位 ,gen(G)
// 下面进行回归
regress 评价量 团购价元 商品毛重kg
// 下面的语句可帮助我们把回归结果保存在Word文档中
// 在使用之前需要运行下面这个代码来安装下这个功能包（运行一次之后就可以注释掉了）
// ssc install reg2docx, all replace
// 如果安装出现connection timed out的错误，可以尝试换成手机热点联网，如果手机热点也不能下载，就不用这个命令吧，可以自己做一个回归结果表，如果觉得麻烦就直接把回归结果截图。
est store m1
reg2docx m1 using m1.docx, replace
// *** p<0.01  ** p<0.05 * p<0.1

// Stata会自动剔除多重共线性的变量
regress 评价量 团购价元 商品毛重kg A1 A2 A3 B1 B2 B3 B4 B5 B6 B7 B8 B9 C1 C2 D1 D2 D3 D4 D5 E1 E2 E3 E4 F1 F2 G1 G2 G3 G4
est store m2
reg2docx m2 using m2.docx, replace

// 得到标准化回归系数
regress 评价量 团购价元 商品毛重kg, b 

// 画出残差图
regress 评价量 团购价元 商品毛重kg A1 A2 A3 B1 B2 B3 B4 B5 B6 B7 B8 B9 C1 C2 D1 D2 D3 D4 D5 E1 E2 E3 E4 F1 F2 G1 G2 G3 G4
rvfplot 
// 残差与拟合值的散点图
graph export a1.png ,replace
// 残差与自变量团购价的散点图
rvpplot  团购价元
graph export a2.png ,replace

// 为什么评价量的拟合值会出现负数？
// 描述性统计并给出分位数对应的数值
summarize 评价量,d

// 作评价量的概率密度估计图
kdensity 评价量 
graph export a3.png ,replace

// 异方差BP检验
estat hettest ,rhs iid

// 异方差怀特检验
estat imtest,white

// 使用OLS + 稳健的标准误
regress 评价量 团购价元 商品毛重kg A1 A2 A3 B1 B2 B3 B4 B5 B6 B7 B8 B9 C1 C2 D1 D2 D3 D4 D5 E1 E2 E3 E4 F1 F2 G1 G2 G3 G4, r
est store m3
reg2docx m3 using m3.docx, replace

// 计算VIF
estat  vif

// 逐步回归（一定要注意完全多重共线性的影响）
// 向前逐步回归（后面的r表示稳健的标准误）
stepwise reg 评价量 团购价元 商品毛重kg A1 A3 B1 B2 B3 B4 B5 B6 B7 B9 C1 D1 D2 D3 D4 E1 E2 E3 F1 G1 G2 G3,  r pe(0.05)
// 向后逐步回归（后面的r表示稳健的标准误）
stepwise reg 评价量 团购价元 商品毛重kg A1 A3 B1 B2 B3 B4 B5 B6 B7 B9 C1 D1 D2 D3 D4 E1 E2 E3 F1 G1 G2 G3,  r pr(0.05)
// 向后逐步回归的同时使用标准化回归系数（在r后面跟上一个b即可）
stepwise reg 评价量 团购价元 商品毛重kg A1 A3 B1 B2 B3 B4 B5 B6 B7 B9 C1 D1 D2 D3 D4 E1 E2 E3 F1 G1 G2 G3,  r b pr(0.05)


// 补充语法 (大家不需要具体的去学Stata软件，掌握我课堂上教给大家的一些命令应对数学建模比赛就可以啦)
// 事实上大家学好Excel，学好后应对90%的数据预处理问题都能解决
// (1) 用已知变量生成新的变量 
generate lny = log(评价量)  
generate price_square = 团购价元 ^2
generate interaction_term = 团购价元*商品毛重kg

// (2) 修改变量名称，因为用中文命名变量名称有时候可能容易出现未知Bug
rename 团购价元 price

lasso回归

LASSO 回归也叫套索回归，是通过生成一个惩罚函数是回归模型中的变量系数进行压缩，达到防止过度拟合，解决严重共线性的问题，LASSO 回归最先由英国人Robert Tibshirani提出，目前在预测模型中应用非常广泛。在新格兰文献中，有大牛提出，对于变量过多而且变量数较少的模型拟合，首先要考虑使用LASSO 惩罚函数。

变量过多会导致多重共线性问题造成回归系数不显著，甚至导致ols估计失效。

image-20240128130519231

Lasso回归的原理

image-20240128132142605

使用lasso回归分析

image-20240128132556845
image-20240128132542976
image-20240128132709587

什么时候使用lasso回归

image-20240128132838334

代码学习与实现：

Stata：拉索开心读懂-Lasso入门 (lianxh.cn)