MATLAB运动学之蒙特卡罗法求积分与机器人工作域分析

蒙特卡罗法又叫做统计模拟法、随机抽样技术,是一种随机模拟方法以概率和统计理论方法为基础的一种计算方法,通俗来说是可以使用随机数来解决很多计算问题的一种方法,很直观简单,尤其对于一些求解积分无解的情况,非常好使且简单粗暴。

蒙特卡罗法求面积(定积分)

y = x² 为例,我们需要求出 x 在[0,10]相对应的 y 在[0,100] 所围成的曲线面积,在我们有了微积分的知识之后,我们可以通过对这个函数的原函数做差来求解(1/3*10³-1/3*0³=1000/3),这种叫做解析解,也就是通过数学公式求出来的解。

除了这种求积分的方法,我们接下来介绍的就是蒙特卡罗法。
将大量随机点散落到整个矩形,然后计算散落在围成曲线下的点的数量的占比就可以得出曲线面积了。
曲线围成的面积=整个矩形区间的面积 * 曲线下方的点的个数的占比

需要注意的是,蒙特卡罗法的前提条件是区间的值要么全是正值,要么全是负值,如果不是的情况就分区再求积分。 

是不是有了这方法,不管什么曲线围成的面积,都不在话下,就这么简单粗暴好用哈哈。

%使用非负整数 seed 为随机数生成函数提供种子,以使 rand、randi 和 randn 生成可预测的数字序列。
rng(0);
set(0,'defaultAxesFontName', 'Monospaced');  % 防止中文乱码
set(gcf, 'position', [10, 20, 1000, 700]);
%f = suptitle('求解y=x^2定积分');
%set(f, 'fontsize', 20); 
L = 10;  % 积分区间长度
fs = 1 / 1e3; % 采样率0.001
x = 0 : fs : L;
y = x .^ 2;  
S = L * (L ^ 2);  %矩形面积,这个示例就是1000% 随机点的数量(作对比)
N_Lis = [10, 100, 1000, 10000];% 解析解(原函数做差值)
res_integ = 1/3 * (10^3 - 0^3); % 近似解
%figure(1); clf;
for n = 1 : length(N_Lis)cnt = 0;x_random = L * rand(1, N_Lis(n));  % 随机点xy_random = L ^ 2 * rand(1, N_Lis(n));  % 随机点y% 统计曲线下面的点的数量for i = 1 : N_Lis(n)if y_random(i) <= x_random(i) ^ 2cnt = cnt + 1;endendres_appro = cnt / N_Lis(n) * S;% 画图对比subplot(2, 2, n);plot(x, y, 'k', 'linewidth', 2); hold on;area(x, y, 'facecolor','c'); hold on;scatter(x_random, y_random, 10, 'r', 'filled', 'markerfacealpha', 0.5);xlabel('x'); ylabel('y'); set(gca, 'fontsize', 14);title(['数学解≈', num2str(res_integ, '%.1f'), '   近似解≈', num2str(res_appro, '%.1f')]);
end

可以看到当随机点从10个增加到10000个的时候,结果对比也可以知道,求出来的这个近似解就越接近解析解(真实值),那么我们在生活当中如果遇到需要求面积的情况,而且连曲线的函数都不清楚的情况下,我们应该知道如何求曲线围成的面积了,比如说,可以撒上一层豆子或者是水,水是最好的(连续,不离散),然后称量下曲线围成的豆子或者水的重量在整个矩形中的占比就可以知道围成的面积了。

无解的情况

有时候求积分是无解的情况,比如下面的三个函数所围成的面积,我们就不能通过数学公式得到解析解或者说非常困难,但是可以快速使用蒙特卡罗法来求其近似解: 

T = 20;
fs = 1 / 1e3;
x0 = -T : fs : T;
y1 = sin(x0.^ 2);
y2 = sin(x0) ./ x0;
y3 = exp(-x0.^2);figure(1); clf;
subplot(3, 1, 1);
plot(x0, y1, 'linewidth', 1.5); ylabel('y'); title('y=sin(x^2)'); set(gca, 'fontsize', 12);
subplot(3, 1, 2);
plot(x0, y2, 'linewidth', 1.5); ylabel('y'); title('y=sin(x)/x'); set(gca, 'fontsize', 12);
subplot(3, 1, 3);
plot(x0, y3, 'linewidth', 1.5); xlabel('x'); ylabel('y'); title('y=e^{-x^2}'); set(gca, 'fontsize', 12);% 绘制围成区域
x = 0 : fs : 2;
y11 = sin(x.^ 2);
y21 = sin(x) ./ x;
y31 = exp(-x.^2);figure(2); clf;
plot(x, y11, 'linewidth', 1.5); hold on;
plot(x, y21, 'linewidth', 1.5); hold on;
plot(x, y31, 'linewidth', 1.5); hold on;
area(x(y11>y31 & y21>y11), y11(y11>y31 & y21>y11), 'facecolor', 'c', 'edgealpha', 0); hold on;
area(x(y11>y31 & y21>y11), y31(y11>y31 & y21>y11), 'facecolor', 'w', 'edgealpha', 0); hold on;
h = legend('y=sin(x^2)', 'y=sin(x)/x', 'y=e^{-x^2}', 'location', 'southwest');
xlabel('x'); ylabel('y'); title('求三条曲线围成的面积'); set(gca, 'fontsize', 12); set(h, 'fontsize', 12);% 蒙特卡罗法求面积
L = 2; 
H = 3;
S = L * H;
N_Lis = [1e1, 1e2, 1e3, 1e4];
figure(3); clf;
for n = 1 : length(N_Lis)N = N_Lis(n);x_random = L * rand(1, N);y_random = H * rand(1, N) - 1;cnt = 0;for i = 1 : Nif (y_random(i) <= sin(x_random(i)^2)) && (y_random(i) <= sin(x_random(i))/x_random(i)) ...&& (y_random(i) >= exp(-x_random(i)^2))cnt = cnt + 1;endendres_appro = cnt / N * S;subplot(2, 2, n);plot(x, y11, 'linewidth', 1.5); hold on;plot(x, y21, 'linewidth', 1.5); hold on;plot(x, y31, 'linewidth', 1.5); hold on;area(x(y11>y31 & y21>y11), y11(y11>y31 & y21>y11), 'facecolor', 'c', 'edgealpha', 0); hold on;area(x(y11>y31 & y21>y11), y31(y11>y31 & y21>y11), 'facecolor', 'w', 'edgealpha', 0); hold on;scatter(x_random, y_random, 10, 'r', 'filled', 'markerfacealpha', 0.5);xlabel('x'); ylabel('y'); title(['样本数=', num2str(N_Lis(n)), '   近似解≈', num2str(res_appro, '%.2f')]); set(gca, 'fontsize', 14); 
endh = suptitle('蒙特卡罗法求图形面积');
set(h, 'fontsize', 18);
set(gcf, 'position', [10, 20, 800, 700]);

只需要将随机点(样本数)增加到基本覆盖整个区域,我们就可以得到所围成的图形里面的样本数的占比,这样就近似求出了这个所围成的面积了。

机器人工作区域

在机器人领域,我们也可以使用蒙特卡罗法模拟出末端执行器的运动区域,这样对于我们关注机器人的所能工作的范围有一个更直观的了解。

%定义D-H参数
a2 = 0.420;
a3 = 0.375;
d2 = 0.138 + 0.024;
d3 =-0.127 -0.024;
d4 = 0.114 + 0.021;
d5 = 0.114 + 0.021;
d6 = 0.090 + 0.021;for i = 1:100000
%角度范围是[-pi,pi],rand返回(0,1) 内均匀分布的随机数
%模拟各关节的角度
theta1 = -pi + 2*pi*rand;
theta2 = 0 + 2*pi*rand;
theta3 =-(5/6)*pi + (5/3)*pi*rand;
theta4 = -pi + 2*pi*rand;
theta5 = -pi + 2*pi*rand;
theta6 = -pi + 2*pi*rand;%XYZ就是关节的末端位置值(不考虑方向)
x(i) = a2*cos(theta1)*cos(theta2)+a3*cos(theta1)*cos(theta2+theta3)-d5*cos(theta1)*sin(theta2+theta3+theta4)-sin(theta1)*(d2+d3+d4)-d6*(cos(theta5)*sin(theta1)-cos(theta1)*cos(theta2+theta3+theta4)*sin(theta5));
01:46y(i) = d6*(cos(theta1)*cos(theta5)+cos(theta2+theta3+theta4)*sin(theta1)*sin(theta5))+a3*sin(theta1)*cos(theta2+theta3)-d5*sin(theta1)*sin(theta2+theta3+theta4)+cos(theta1)*(d2+d3+d4)+a2*cos(theta2)*sin(theta1);z(i) =-a3*sin(theta2+theta3)-a2*sin(theta2)-d5*cos(theta2+theta3+theta4)-d6*sin(theta5)*sin(theta2+theta3+theta4);
endplot3(x,y,z,'b.','MarkerSize',0.5)

我们这里让机器人的关节随机运行10万次,也就是10万个随机点,通过plot3函数,画出这个六轴机械臂末端执行器所处空间的能够工作的范围了,基本上可以看到能够覆盖机器人所能够工作的区域了。

也可以观察XY组成的侧面,或者另外两根轴组成的侧面情况

其中中间白色圆心部分,是机械臂末端所不能运动到的地方。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/216105.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

书-二分查找找某个数字p155

#include<stdio.h> int main(){int a[10]{1,4,5,6,7,8,23,34,90,14567};int mid;int low0;int high9;while(low<high){mid(lowhigh)/2;//数组分成两段&#xff0c;前一段low-mid&#xff0c;后一段mid-highif (a[mid]<23)//因为已经是排序好的了&#xff0c;所以如…

基于ssm志愿者招募网站源码和论文

网络的广泛应用给生活带来了十分的便利。所以把志愿者招募管理与现在网络相结合&#xff0c;利用java技术建设志愿者招募网站&#xff0c;实现志愿者招募的信息化。对于进一步提高志愿者招募管理发展&#xff0c;丰富志愿者招募管理经验能起到不少的促进作用。 志愿者招募网站…

虹科分享 | CanEasy多场景应用,让汽车总线测试更简单

CanEasy是一个基于Windows的总线工具&#xff0c;用于分析和测试CAN、CAN FD和LIN以及汽车以太网系统。通过高度自动化和简单的配置模拟总线流量&#xff0c;CanEasy可用于分析真实网络、模拟虚拟系统&#xff0c;以及在整个开发过程中进行剩余总线模拟&#xff0c;实现从测试到…

【力扣100】48.旋转图像

添加链接描述 class Solution:def rotate(self, matrix: List[List[int]]) -> None:n len(matrix)# Python 这里不能 matrix_new matrix 或 matrix_new matrix[:] 因为是引用拷贝matrix_new [[0] * n for _ in range(n)]for i in range(n):for j in range(n):matrix_ne…

“2024杭州国际安防产品展览会”加快推进浙江平安城市体系现代化

2024杭州国际安防产品展览会&#xff0c;将于2024年3月份在杭州国际博览中心隆重召开。作为全球安防产业的重要盛会&#xff0c;该展览会以其“专业化、市场化”为特点&#xff0c;聚焦每年核心主题&#xff0c;以专业化为核心&#xff0c;打造品牌项目&#xff0c;使之成为备受…

Ansys Lumerical | 采用一维光栅的出瞳扩展器的优化

附件下载 联系工作人员获取附件 本文演示了一种仿真方法&#xff0c;并举例说明了使用一维光栅的出瞳扩张器&#xff08;EPE&#xff09;系统的优化示例。 在此工作流程中&#xff0c;我们使用 Lumerical 构建光栅模型&#xff0c;并使用 RCWA 求解器模拟其响应。完整的EPE系…

一元四次方程解法,亲测ok

一元四次方程求解方法及程序 - 知乎 (zhihu.com) 自己写的不行&#xff0c;这个博主的完全ok&#xff0c;赞&#xff01;&#xff01;&#xff01;

关于上传的图片-File,Base64,Blob之间的转换

base64转file export function base64toFile (dataurl, filename) {let arr dataurl.split(,)let mime arr[0].match(/:(.*?);/)[1]let suffix mime.split(/)[1]let bstr atob(arr[1])let n bstr.lengthlet u8arr new Uint8Array(n)while (n--) {u8arr[n] bstr.charC…

arcgis api for js 中使用API的代理页面(跨越配置)

以下仅作为自己阅读官网api的对reques的理解做的备忘笔记。一知半解&#xff0c;仅供参考。 1、获取或者构建第三方代理 官网解释&#xff1a;代理在其自己的 Web 服务器上安装并运行&#xff0c;而不是在 Esri 服务器或安装了 ArcGIS Enterprise 的计算机上安装和运行&#…

(新手必看)自定义数据传输通信协议+STM32代码详解

前言 本篇博客主要学习和了解一些单片机协议的格式&#xff0c;在对传输大数据或者要求准确性的时候&#xff0c;都需要通过协议来发送接收&#xff0c;下面通过了解协议的基本构成和代码来分析和实现协议的发送和接收。本篇博客大部分是自己收集和整理&#xff0c;如有侵权请联…

机器学习-KL散度的直观理解+代码

KL散度 直观理解&#xff1a;KL散度是一种衡量两个分布之间匹配程度的方法。通常在概率和统计中&#xff0c;我们会用更简单的近似分布来代替观察到的数据或复杂的分布&#xff0c;KL散度帮我们衡量在选择近似值时损失了多少信息。 在信息论或概率论中&#xff0c;KL散度&#…

Java+Swing: 从数据库中查询数据并显示在表格中 整理11

分析&#xff1a;要想从数据库中查询数据并分页展示到表格中&#xff0c;我觉得应该按照这个思路&#xff1a;首先就是发起请求&#xff0c;此时需要向数据库中传递三个参数&#xff1a;当前页码&#xff08;pageNum&#xff09;、每一页的数量&#xff08;pageSize&#xff09…

好用的API接口大汇总,含免费次数

AI绘画-Stable Diffusion&#xff1a;通过AI 生成图片&#xff0c;包括图生文、文生图等。AI绘画-Mid Journey&#xff1a;使用 Midjourney 目前全球领先的图片大模型&#xff0c;其能根据输入文字提供极其优秀的AI绘画作品。运营商三要素 API&#xff1a;输入姓名、身份证号码…

Java八股文面试全套真题【含答案】- AJAX Axios篇

AJAX 是什么&#xff1f;它的全称是什么&#xff1f; 答案&#xff1a;AJAX 是 Asynchronous JavaScript and XML&#xff08;异步 JavaScript 和 XML&#xff09;的缩写。它是一种用于在后台与服务器进行数据交换的技术&#xff0c;实现异步加载数据而无需刷新整个页面。AJAX …

TypeScript 第十一节:命名空间

一、命名空间 TypeScript 中的命名空间&#xff08;Namespace&#xff09;用于将代码组织到逻辑分组中。在 TypeScript 中&#xff0c;命名空间是一个独立作用域中的代码集合。 1、示例 下面是一个简单的 TypeScript 命名空间示例&#xff1a; namespace MyNamespace {export…

Knowledge Graph知识图谱—8. Web Ontology Language (OWL)

8. Web Ontology Language (OWL) 在RDFs不可能实现&#xff1a; Property cardinalities, Functional properties, Class disjointness, we cannot produce contradictions, circumvent the Non Unique Naming Assumption, circumvent the Open World Assumption 8.1 OWL Tr…

window下make无法识别,同时缺少mingw32-make.exe文件

解决方法&#xff1a;window下make无法识别&#xff0c;同时缺少mingw32-make.exe文件_无法将make项识别为cmdlet-CSDN博客

网工排查网络故障,有这两款软件就够了

网络工程师的工作中&#xff0c;排查网络故障占很大一部分。领导让你查网络故障&#xff0c;批量Ping几十台电脑&#xff0c;结果你Ping了一个小时还没好。你总不能说批量Ping就是这么慢吧&#xff1f; 所以&#xff0c;在这里介绍两款网络工程师常用的排查网络故障工具。 Qu…

K8S学习指南(2)-docker的基本使用

文章目录 引言安装 DockerDocker 基本概念1. 镜像&#xff08;Images&#xff09;示例&#xff1a;拉取并运行一个 Nginx 镜像 2. 容器&#xff08;Containers&#xff09;示例&#xff1a;查看运行中的容器 3. 仓库&#xff08;Repository&#xff09;示例&#xff1a;推送镜像…

知名火锅连锁企业,IT 团队如何在数千家门店中先于用户发现故障

该知名火锅连锁企业是中国领先的餐饮企业&#xff0c;上千家门店遍布全球&#xff0c;由于门店餐饮行业的特殊性&#xff0c;需要靠前部署服务&#xff0c;所以在每家餐厅中&#xff0c;会部署相应的服务器&#xff0c;及相应 IT 设备&#xff0c;本地会运行POS、会员、下单等业…