非线性规划规划问题
线性规划问题的复习
在学习非线性规划问题的过程中学到了蒙特卡洛模拟。
对于蒙特卡洛模拟的思路就是：
1：从约束条件中求出每个变量的范围
2：在上述范围中用随机生成若干组实验点，并验证他们是否满足所有的约束条件，若满足，则将其划分到可行组，之后从可行组中找到函数的最大值或最小值
注意：每一个决策变量都是有界的，才能用蒙特卡洛模拟出来随机数。
根据约束条件有时可以减少决策变量的个数
放朔技巧
由上述计算过程可以得到x1和x3的可能取值的范围，分别是[20, 30], [-10, 16]。由此结果我们可以生成很多随机数。思想是：只要我生成的随机数足够多酒呢能找到最大值。但是由于耗时的原因。模拟的次数要把握好。
% 预备知识
% (1) format long g 可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
5/7
5895*514100
format long g
5/7
5895*514100
% (2)unifrnd(a,b,m,n)可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。(等价于 a + rand(m,n)*(b-a))
unifrnd(0,5,4,3)
% 4.07361843196589 3.16179623112705 4.78753417717149
% 4.5289596853781 0.487702024997048 4.82444267599638
% 0.63493408146753 1.39249109433524 0.788065408387741
% 4.5668792806951 2.73440759602492 4.85296390880308
% 正式代码
clc,clear;
tic %开始计时
n=10000000; %生成的随机数组数
x1=unifrnd(20,30,n,1); % 生成在[20,30]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(-10,16,n,1); % 生成在[-10,16]之间均匀分布的随机数组成的n行1列的向量构成x3
fmax=-inf;
for i=1:n
x = [x1(i), x2(i), x3(i)]; %构造x向量
if (-x(1)+2*x(2)+2*x(3)>=0) & (x(1)+2*x(2)+2*x(3)<=72) % 约束条件
result = x(1)*x(2)*x(3); % 如果满足条件就计算函数值
if result > fmax
fmax = result;
X = x; % 并且将此时的x1 x2 x3保存到一个变量中
end
end
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计时结束
缩小范围重新模拟得到更加精确的取值
clc,clear;
tic
n=10000000; %生成的随机数组数
x1=unifrnd(22,23,n,1); %由于第一次模拟生成的结果缩小x1和x3的范围
x3=unifrnd(11,13,n,1);
fmax=-inf;
for i=1:n
x = [x1(i), x2(i), x3(i)];
if (-x(1)+2*x(2)+2*x(3)>=0) & (x(1)+2*x(2)+2*x(3)<=72)
result = x(1)*x(2)*x(3);
if result > fmax
fmax = result;
X = x;
end
end
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc
书店买书问题（0-1规划）
当有n家店，m本书，那么所有可能的情况有
nm种，所有当nm非常大的时候我们再用蒙特卡洛模拟就不现实，要使用其他的方法。
% 预备知识
% (1)unique函数: 剔除一个矩阵或者向量的重复值，并将结果按照从小到大的顺序排列
% adj. 唯一的; 独一无二的 [ju'ni:k]
unique([1 2 5; 6 8 9;2 4 6])
% ans = [1;2;3;4;5;6;7;8;9]
unique([5 6 8 8 4 1 6 2 2 4 8 4 5 6])
% ans = [1 2 4 5 6 8]
% (2)randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
randi([-5,5],2,6)
min_money = +Inf;
min_result = randi([1, 6],1,5);
n = 100000; % 蒙特卡罗模拟的次数
M = [18 39 29 48 59
24 45 23 54 44
22 45 23 53 53
28 47 17 57 47
24 42 24 47 59
27 48 20 55 53]; %售价
freight = [10 15 15 10 10 15]; % 第i家店的运费
for k = 1:n
result = randi([1, 6],1,5); % 在1-6这些整数中随机抽取一个1*5的向量，表示这五本书分别在哪家书店购买
index = unique(result); % 在哪些商店购买了商品，因为我们等下要计算运费
money = sum(freight(index)); % 计算买书花费的运费
% 计算总花费：刚刚计算出来的运费 + 五本书的售价
for i = 1:5
money = money + M(result(i),i);
end
if money < min_money
min_money = money
min_result = result
end
end
min_money % 18+39+48+17+47+20
min_result
导弹追踪问题
因为导弹始终对准B点，所以在任意的时刻导弹的行进轨迹的切线都要经过B船的位置
将连续的模型离散化：将时间T划分为多个▲T
% 预备知识
% mod(m, n)表示m/n的余数
mod (8,3)
% ans = 2
% 设置横纵坐标的范围并标上字符
x = 1:0.01:3;
y = x .^ 2;
plot(x,y) % 画出x和y的图形
axis([0 3 0 10]) % 设置横坐标范围为[0, 3] 纵坐标范围为[0, 10]
pause(3) % 暂停3秒后再继续接下来的命令
text(2,4,'清风') % 在坐标为(2,4)的点上标上字符串：清风
close % 关闭图形窗口
只是对结果进行的计算
其核心是先计算出一个▲T之后的B船的位置，之后由B船的位置可以得到导弹的实时方向，再通过其方向计算出导弹的位置
之后每一次经过一个▲T就判断导弹射程是否已经满了，导弹与B船的位置是否足够小，来决定是否结束循环。
clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2; % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
while(dd>=0.001) % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）
t=t+dt; % 更新导弹击落B船的时间
d=d+3*v*dt; % 更新导弹飞行的距离
x(2)=20+t*v*m; y(2)=t*v*m; % 计算新的B船的位置 （注：m=sqrt(2)/2）
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 更新导弹与B船的距离
tan_alpha=(y(2)-y(1))/(x(2)-x(1)); % 计算斜率，即tan(α)
cos_alpha=sqrt(1/(1+tan_alpha^2)); % sec(α)^2 = (1+tan(α)^2)
sin_alpha=sqrt(1-cos_alpha^2); % sin(α)^2 +cos(α)^2 = 1
x(1)=x(1)+3*v*dt*cos_alpha; y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置
if d>50 % 导弹的有效射程为50个单位
disp('导弹没有击中B船');
break; % 退出循环
end
if d<=50 & dd<0.001 % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）
disp(['导弹飞行',num2str(d),'单位后击中B船'])
disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
end
end
在导弹追击的过程中画出示意图。
clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);
for i=1:2
plot(x(i),y(i),'.k','MarkerSize',1); % 画出导弹和B船所在的坐标，点的大小为1，颜色为黑色(k)，用小点表示
grid on;
hold on;
axis([0 30 0 10])
k = 0; % 引入一个变量 为了控制画图的速度（因为Matlab中画图的速度超级慢）
while(dd>=0.001) % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）
t=t+dt; % 更新导弹击落B船的时间
d=d+3*v*dt; % 更新导弹飞行的距离
x(2)=20+t*v*m; y(2)=t*v*m; % 计算新的B船的位置 （注：m=sqrt(2)/2）
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 更新导弹与B船的距离
tan_alpha=(y(2)-y(1))/(x(2)-x(1)); % 计算斜率，即tan(α)
cos_alpha=sqrt(1/(1+tan_alpha^2)); % 利用公式：sec(α)^2 = (1+tan(α)^2) 计算出cos(α)
sin_alpha=sqrt(1-cos_alpha^2); % 利用公式： sin(α)^2 +cos(α)^2 = 1 计算出sin(α)
x(1)=x(1)+3*v*dt*cos_alpha; y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置
k = k +1 ;
if mod(k,500) == 0 % 每刷新500次时间就画出下一个导弹和B船所在的坐标
for i=1:2
plot(x(i),y(i),'.k','MarkerSize',1);
hold on; % 不关闭图形，继续画图
end
pause(0.001); % 暂停0.001s后再继续下面的操作
end
if d>50 % 导弹的有效射程为50个单位
disp('导弹没有击中B船');
break; % 退出循环
end
if d<=50 & dd<0.001 % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）
disp(['导弹飞行',num2str(d),'个单位后击中B船'])
disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
end
end
旅行商问题（TSP）
由数学归纳法容易知到当有n个城市的时候，有（n-1）！种路径可以走。
由于阶乘的增长程度十分大，所以当城市数目超过10的时候用此种方法就不现实。要换另一种方式来做。
% 预备知识
% randperm函数的用法
randperm(5) % 生成1-5组成的一个随机序列(类似于洗牌的操作)
% 3 5 1 2 4
% 1 4 5 3 2
prod函数是连续求乘积的函数
prod(1:5)
% ans = 120
clear;clc
% 只有10个城市的简单情况
coord =[0.6683 0.6195 0.4 0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488;0.2536 0.2634 0.4439 0.1463 0.2293 0.761 0.9414 0.6536 0.5219 0.3609]';
n = size(coord,1); % 城市的数目
figure(1) % 新建一个编号为1的图形窗口
plot(coord(:,1),coord(:,2),'o'); % 画出城市的分布散点图
for i = 1:n
text(coord(i,1)+0.01,coord(i,2)+0.01,num2str(i)) % 在图上标上城市的编号（加上0.01表示把文字的标记往右上方偏移一点）
end
hold on % 等一下要接着在这个图形上画图的
d = zeros(n); % 初始化两个城市的距离矩阵全为0
for i = 2:n
for j = 1:i
coord_i = coord(i,:); x_i = coord_i(1); y_i = coord_i(2); % 城市i的横坐标为x_i，纵坐标为y_i
coord_j = coord(j,:); x_j = coord_j(1); y_j = coord_j(2); % 城市j的横坐标为x_j，纵坐标为y_j
d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2); % 计算城市i和j的距离
end
end
d = d+d'; % 生成距离矩阵的对称的一面
min_result = +inf; % 假设最短的距离为min_result，初始化为无穷大，后面只要找到比它小的就对其更新
min_path = [1:n]; % 初始化最短的路径就是1-2-3-...-n
N = 10000000; % 蒙特卡罗模拟的次数
for i = 1:N % 开始循环
result = 0; % 初始化走过的路程为0
path = randperm(n); % 生成一个1-n的随机打乱的序列
for i = 1:n-1
result = d(path(i),path(i+1)) + result; % 按照这个序列不断的更新走过的路程这个值
end
result = d(path(1),path(n)) + result; % 别忘了加上从最后一个城市返回到最开始那个城市的距离
if result < min_result % 判断这次模拟走过的距离是否小于最短的距离，如果小于就更新最短距离和最短的路径
min_path = path;
min_result = result
end
end
min_path
min_path = [min_path,min_path(1)]; % 在最短路径的最后面加上一个元素，即第一个点（我们要生成一个封闭的图形）
n = n+1; % 城市的个数加一个（紧随着上一步）
for i = 1:n-1
j = i+1;
coord_i = coord(min_path(i),:); x_i = coord_i(1); y_i = coord_i(2);
coord_j = coord(min_path(j),:); x_j = coord_j(1); y_j = coord_j(2);
plot([x_i,x_j],[y_i,y_j],'-') % 每两个点就作出一条线段，直到所有的城市都走完
pause(0.5) % 暂停0.5s再画下一条线段
hold on
end
总结
经过这些个样例，显而易见的，蒙特卡洛模拟就是暴力，只要我的数据够大，就没有我解不出来的模型，较为简单。但是我们在应用中要结合实际问题的复杂程度来确定需要蒙特卡洛模拟的次数，作答。
扫一扫
970
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
