最近发现符号计算库sympy的强大之处在于公式推导,然后帮小伙伴解决了一个线性方程组求解的问题,特此记录一下。
另外,表扬csdn支持一键导入markdown,这样就可以将从Anaconda导出的md一键导入到日志编辑框,非常省事!
问题描述:
import numpy as np
import sympy as sp
import scipy.optimize as opt
# If all you want is the best pretty printing, use the init_printing() function.
# This will automatically enable the best printer available in your environment.
# 自动使用最漂亮的输出方式
sp.init_printing()
sumpy.init_printing() 自动使用最漂亮的公式输出方式
1. 利用 sympy 建立方程组,推导公式
2. 将 sympy 推导出来的公式转化为代码
3.利用 scipy.optimize.least_squares()函数求方程的根
1. 利用 sympy 建立方程组,推导公式
1.1 构造矩阵
# unknown parameters; to be solved
theta1 = sp.Symbol('theta1', real=True)
theta2 = sp.Symbol('theta2', real=True)
phi = sp.Symbol('phi', real=True)
# known hyperparameters; to be set
psi1 = sp.Symbol('psi1', real=True)
psi2 = sp.Symbol('psi2', real=True)
alpha = sp.cos(psi1) # real number
beta = sp.sin(psi1) * sp.exp(psi2 * 1j) # complex number
# construct the matrices
x = theta1 * 2
y = theta2 * 2
# sumpy的 I 等于 python自带的 1j
A = sp.Matrix([[1 + 1j * sp.cos(x), 1j * sp.sin(x)], [1j * sp.sin(x), 1 - 1j * sp.cos(x)]]) * np.sqrt(1/2)
B = sp.Matrix([[sp.cos(y), sp.sin(y)],[sp.sin(y), - sp.cos(y)]]) * 1j
C = sp.Matrix([[1], [0]])
all_phase = sp.exp( phi*1j )
D = A * B * C * all_phase
D1 = sp.simplify(D)
D1
alpha,beta
1.2 定义方程并化简
sumpy.simplify(expr) 可以将表达式 expr 化到最简
J1 = sp.simplify(sp.re(D1[0]) - sp.re(alpha))
J2 = sp.simplify(sp.im(D1[0]) - sp.im(alpha))
J3 = sp.simplify(sp.re(D1[1]) - sp.re(beta))
J4 = sp.simplify(sp.im(D1[1]) - sp.im(beta))
分别看4个方程
J1
J2
J3
J4
1.3 求雅可比行列式
因为待求参数只有3个,因此只需要联立3个方程。
因为 alpha为实数,这里选择J1、J3和J4联立方程组,并分别求出它们对各个参数的偏导数,用于构造雅可比行列式。
Jacobian 雅可比行列式,以n个n元函数的偏导数为元素的行列式。
transcendental equation 超越方程
dJ1_theta1 = sp.simplify(sp.diff(J1, theta1))
dJ1_theta2 = sp.simplify(sp.diff(J1, theta2))
dJ1_phi = sp.simplify(sp.diff(J1, phi))
dJ3_theta1 = sp.simplify(sp.diff(J3, theta1))
dJ3_theta2 = sp.simplify(sp.diff(J3, theta2))
dJ3_phi = sp.simplify(sp.diff(J3, phi))
dJ4_theta1 = sp.simplify(sp.diff(J4, theta1))
dJ4_theta2 = sp.simplify(sp.diff(J4, theta2))
dJ4_phi = sp.simplify(sp.diff(J4, phi))
雅可比行列式的每一行
dJ1_theta1,dJ1_theta2,dJ1_phi
dJ3_theta1,dJ3_theta2,dJ3_phi
dJ4_theta1,dJ4_theta2,dJ4_phi
2. 将 sympy 推导出来的公式转化为代码
x=[theta1, theta2, phi];
f consists of [J1,J3,J4]
为方程组, 为
# x:[theta1, theta2, phi]; psi=pi/4
# f consists of [J1,J3,J4]
# 目标函数
def fun_tf_ls(x, psi1, psi2):
f =[ - np.sqrt(1/2) * np.sin(x[2]) * np.cos(2 * x[1])
- np.sqrt(1/2) * np.cos(x[2]) * np.cos(2 * (x[0] - x[1])) - np.cos(psi1),
- np.sqrt(1/2) * np.sin(x[2]) * np.sin(2 * x[1])
- np.sqrt(1/2) * np.sin(2 * (x[0] - x[1])) * np.cos(x[2]) - np.sin(psi1) * np.cos(psi2),
- np.sqrt(1/2) * np.sin(x[2]) * np.sin(2 * (x[0] - x[1]))
+ np.sqrt(1/2) * np.sin(2 * x[1]) * np.cos(x[2]) - np.sin(psi1) * np.sin(psi2) ] # 3个方程
return f
# 雅可比行列式
def deri_tf_ls(x, psi1, psi2):
df = np.array([[2 * np.sqrt(1/2) * np.sin(2 * (x[0] - x[1])) * np.cos(x[2]),
2 * np.sqrt(1/2) * ( np.sin(x[2])*np.sin(2*x[1])-np.sin(2*(x[0]-x[1]))*np.cos(x[2])),
np.sqrt(1/2) * (np.sin(x[2]) * np.cos(2*(x[0]-x[1]))-np.cos(x[2])*np.cos(2 * x[1]))],
[-2 * np.sqrt(1/2) * np.cos(x[2]) * np.cos(2 * (x[0] - x[1])),
-2 * np.sqrt(1/2) * (np.sin(x[2]) * np.cos(2*x[1])-np.cos(x[2])*np.cos(2*(x[0]-x[1]))),
np.sqrt(1/2)*(np.sin(x[2]) * np.sin(2 *(x[0]-x[1]))-np.sin( 2 * x[1])*np.cos(x[2]))],
[-2 * np.sqrt(1/2) * np.sin(x[2]) * np.cos(2 * (x[0] - x[1])),
2 * np.sqrt(1/2) * (np.sin(x[2])*np.cos(2*(x[0]-x[1])) + np.cos(x[2])*np.cos(2*x[1])),
-np.sqrt(1/2) * (np.sin(x[2])*np.sin(2*x[1]) + np.cos(x[2])*np.sin(2*(x[0]-x[1])))]
]) # 3 X 3
return df
3.利用 scipy.optimize.least_squares()函数求方程的根
scipy.optimize.least_squares(fun, x0, jac=‘2-point’, bounds=(-inf, inf), method=‘trf’, ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, loss=‘linear’, f_scale=1.0, diff_step=None, tr_solver=None, tr_options={}, jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs={})
重要参数:method 和 loss
method指定优化算法。对于有界优化问题,trf(Trust Region Reflective algorithm)是最健壮的。
loss 可以指定代价的计算形式,'linear’指标准的最小二乘代价。
opt.least_squares()建议先算出雅可比行列式的解析表达式,然后再传给jac参数,否则函数会使用有限差分算法估计雅可比行列式,这样计算精度比较低,而且用时更长。
psi1_0 = np.pi/4
psi2_0 = np.pi/2
x0 = np.array([1,1,1])
# 限定[theta1, theta2, phi]的定义域为[0,pi]
# bounds=(lower_bound, upper_bound);
# lower_bound和upper_bound可以为具体数值,也可以为np.inf(正无穷或-np.inf(负无穷)
# 给每个自变量单独指定定义域:bounds=([0,0,0], [np.pi, np.pi, np.pi])
# 为所有自变量指定相同的定义域: bounds=(0,np.pi)
sol_tf = opt.least_squares(fun_tf_ls, x0, args=(psi1_0,psi2_0), jac=deri_tf_ls, bounds=(0, np.pi))
sol_tf
active_mask: array([0, 0, 0])
cost: 2.56332266214165e-13
fun: array([-6.88199085e-07, -6.07527060e-10, -1.97601070e-07])
grad: array([ 8.68936726e-10, -2.15471831e-09, -3.40677099e-10])
jac: array([[-8.69567557e-04, 1.73913633e-03, 1.11924532e-03],
[-8.41719879e-07, 1.41421387e+00, -7.07106584e-01],
[-1.36892283e-03, 4.99355034e-04, -6.07527016e-10]])
message: '`gtol` termination condition is satisfied.'
nfev: 11
njev: 11
optimality: 3.385286057228377e-09
status: 1
success: True
x: array([0.78557471, 1.57048889, 1.57018145])
3.1 方程的根
parameters_final = sol_tf.x
parameters_final
array([0.78557471, 1.57048889, 1.57018145])
3.1.1 以 pi 为单位
parameters_final_pi = sol_tf.x / np.pi
parameters_final_pi
array([0.2500562 , 0.49990214, 0.49980428])
3.1.2 以 ° 为单位
parameters_final_degree = sol_tf.x / np.pi * 180
parameters_final_degree
array([45.01011548, 89.98238504, 89.96477012])
3.2 目标函数的最终值(尽可能接近0)
3.2.1 scipy.optimize.least_squares()的结果
sol_tf.fun
array([-6.88199085e-07, -6.07527060e-10, -1.97601070e-07])
3.2.2 根据sympy模型正向计算的结果
psi = pi/4, [theta1,theta2,phi] = [pi/4, pi/2, pi/2]
J_total = sp.Matrix([J1,J3,J4])
J_total_test = J_total.evalf(subs={theta1:parameters_final[0], theta2:parameters_final[1], phi:parameters_final[2],
psi1:psi1_0, psi2:psi2_0}, chop=True)
J_total_test
忽略numpy、scipy和sympy的计算精度差异,两者得到的结果是一致的,最终目标函数都被优化到0,说明方程组得解。