pid算法c语言详解：计算PIπ值
pid算法c语言详解：计算PIπ值可以看到，x、y的值有时趋势完全一样，π的估值则明显有周期性，大概一分钟。了解博途随机数生成原理的应该知道，它是通过读取时间里的纳秒经过换算得到的。π的估值#n1 := 0; // 程序变量初始化 FOR #i := 0 TO #n DO // 利用LGF库内的随机实数程序，生成在0.0-1.0范围内的实数并赋值 #x := "LGF_RandomRange_Real"(minValue := 0.0 maxValue := 1.0 error => "Tag_17" status => "Tag_18" subfunctionStatus => "Tag_12"); #y := "LGF_RandomRange_Real"(minValue := 0.0 max
计算π的方法有很多，前两天看到一个有意思的方法：蒙特卡洛方法（Monte Carlo method）。其实这种方法由来已久，它的思想是根据某项问题建立一个概率事件，通过随机采样来对实际值进行估算。
它的一个简单示例就是估算PI值。如下图：构造一个正方形和一个1/4圆，在整个区域内随机投入点，根据投入点到原点的距离判断是落在圆内还是在圆外。落在圆内的概率其实就是1/4圆与正方形面积比：π/4。
博途仿真结果
样本数量n为100000
程序实现
新建FC块，添加变量
#n1 := 0; // 程序变量初始化
FOR #i := 0 TO #n DO
// 利用LGF库内的随机实数程序，生成在0.0-1.0范围内的实数并赋值
#x := "LGF_RandomRange_Real"(minValue := 0.0 maxValue := 1.0 error => "Tag_17" status => "Tag_18" subfunctionStatus => "Tag_12");
#y := "LGF_RandomRange_Real"(minValue := 0.0 maxValue := 1.0 error => "Tag_19" status => "Tag_13" subfunctionStatus => "Tag_20");
// 样本如果处于圆内，则计数加1
IF (#x*#x #y*#y)<1.0 THEN
#n1 = 1;
END_IF;
END_FOR;
#Pi := 4.0*DINT_TO_LREAL(#n1)/DINT_TO_LREAL(#n); //π=4*n1/n
投入点坐标x、y由伪随机实数生成程序"LGF_RandomRange_Real"生成。官方提供了包含该程序的LGF库(Library of General Functions)。百度网盘：https://pan.baidu.com/s/1rqweR87IHYpPkIRzi_MTaQ 提取码：1210
仿真曲线：
使用博途仿真的结果其实并不稳定，以下分别坐标x、y和π的曲线：
坐标x、y的取值
π的估值
可以看到，x、y的值有时趋势完全一样，π的估值则明显有周期性，大概一分钟。了解博途随机数生成原理的应该知道，它是通过读取时间里的纳秒经过换算得到的。
这其实也验证了：PLC乃至计算机不会产生绝对随机的随机数，只能产生“伪随机数”，这里的“伪”是有规律的意思。
还有一点，同样的程序，在不同的电脑上π的仿真结果并不一致（公司：π在3.22附近；家里：在3.04附近）。搜索了一下相关问题，有这么句话"仿真是基于操作系统的，它的循环时间与实际PLC硬件运行时并不一样。" 难道是因为电脑配置的不同？受循环周期的影响？
使用Python：
同样的程序转到Python上运行一下，误差稳定且小了很多。猜测可能随机数的选取或者程序运行流程更加合理吧。
有兴趣的可以下载个软件试试，这里使用的是PyCharm，纯计算，没啥复杂度，我也属于现学现用。
from random import random
from math import sqrt
n = 100000 # 样本总量
n1 = 0 # 落入圆内数量
for i in range(1 n):
x = random();
y = random();
if sqrt(x**2 y**2) <= 1.0:
n1 = n1 1
pi = 4.0 * (n1/n)
print('Pi的值为 %s' % pi)
样本量同样是100000
总结与后续：
以上的蒙特卡洛方法算是一种统计模拟方法，在求解微分方程，多重积分，特征值等方面应用也很多。它将确定性问题转化成了随机性问题，采样越多，越近似最优解，但永远不是最优解，某种程度上丧失了精确性。
以下是一些关于π的公式（使用PY的计算结果就不贴出来了）：
莱布尼茨公式：
arctan(1) = π/4，带入1即可
sum=0
i = 0
for i in range(0 1000000):
if i % 2 == 0:
sum = 1/(2*i 1)
else:
sum -= 1/(2*i 1)
pi=4*sum
print('Pi的值为 %s' % pi)Machin_梅钦(马青)公式：
import math
pi=16*math.atan(1/5)-4*math.atan(1/239)
print('Pi的值为 %s' % pi)
直接使用现成的反三角函数
sum=0
i = 0
b = 1
for i in range(0 5):
sum = 16*(1/5)**(2*i 1)/(2*i 1)*b-4*(1/239)**(2*i 1)/(2*i 1)*b
b=-b
pi=sum
print('Pi的值为 %s' % pi)欧拉：
sum=0
i = 0
for i in range(1 1000000):
sum = 1/(i*i)
pi=(6*sum)**0.5
print('Pi的值为 %s' % pi)拉玛努金：
k=0时，π≈9801/1103/2/2**0.5=3.1415927300133055；
k=1时，π≈9801/2/2**0.5/(1103 24*(1103 26390)/396!)=3.1415926524195665；
这个我们欣赏一下就行了。
改进型：
k=0时，π≈10005**0.5*426880/13591409=3.141592653589734；
............
高斯-勒让德-迭代算法
迭代的过程对于计算机来说简直太适合了。
import math
a=1
b=1/(2**0.5)
t=0.25
p=1
i=0
for i in range(0 3):
a1 = (a b)/2
b1 = (a*b)**0.5
t1 = t-p*((a-(a b)/2)**2)
p1 = 2*p
pi = ((a1 b1)**2)/(4*t1)
a = a1
b = b1
t = t1
p = p1
i=i 1
print(f'迭代{i}次的π值为：{pi}')
迭代1次的值为：3.1405792505221686
迭代2次的值为：3.141592646213543
迭代3次的值为：3.141592653589794神奇的π：
欧拉恒等式：
爱因斯坦场方程：
黑洞温度（霍金）：
人生苦短，当有所钟！
