我正在研究一个Python优化问题，我需要强制执行一个约束，即两个向量（由点JG和LA形成）之间的角度必须小于160度。在优化过程中应遵守此约束，特别是对于点的初始位置。然而，尽管实施了基于惩罚的方法并使用具有显式非线性约束的 SLSQP 优化器，但角度并未得到考虑。我正在寻求有关如何改进代码的建议，以确保在整个优化过程中严格执行此约束。
我定义了一个惩罚函数，如果在初始位置优化期间角度超过 160 度，则为目标函数添加较大的惩罚。 我还尝试将 NonlinearConstraint 与 SLSQP 优化器结合使用来显式强制执行约束。 尽管做出了这些努力，优化结果有时仍然违反角度约束，特别是在初始配置中。
我希望优化过程能够严格执行角度约束，确保 JG 和 LA 之间的角度保持在 160 度以下，特别是对于初始位置。
在确定点的初始位置后，根据 BL 和 GJ 的最佳角度以及距离的初始猜测（dist_BF、dist_FG、dist_GJ 等）计算距离 dist_AI。这是计算 dist_AI 的代码的特定部分：
points_initial = compute_points(angle_BL_opt, angle_GJ_opt, *initial_guess) dist_AI = np.sqrt((points_initial['A'][0] - points_initial['I'][0])**2 + (points_initial['A'][1] - points_initial['I'][1])**2)
代码：
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve, minimize, minimize_scalar, NonlinearConstraint
# Helper functions
def deg_to_rad(angle_deg):
return angle_deg * np.pi / 180
def rad_to_deg(angle_rad):
return angle_rad * np.pi / 180
def calculate_distance_CE(angle_BL_deg, dist_CE_desired):
angle_BL = deg_to_rad(angle_BL_deg)
points = {
'A': (0, 0),
'B': (-60, -255),
'C': (-4, -669.5)
}
dist_BL = 2110
points['L'] = (points['B'][0] + dist_BL * np.cos(angle_BL), points['B'][1] + dist_BL * np.sin(angle_BL))
dist_BD = 1179.5
points['D'] = (points['B'][0] + dist_BD * np.cos(angle_BL), points['B'][1] + dist_BD * np.sin(angle_BL))
dist_BF = 1663
points['F'] = (points['B'][0] + dist_BF * np.cos(angle_BL), points['B'][1] + dist_BF * np.sin(angle_BL))
angle_DE = angle_BL + np.pi / 2
dist_DE = 30
points['E'] = (points['D'][0] + dist_DE * np.cos(angle_DE), points['D'][1] + dist_DE * np.sin(angle_DE))
if points['E'][1] <= points['D'][1]:
points['E'] = (points['D'][0] - dist_DE * np.cos(angle_DE), points['D'][1] - dist_DE * np.sin(angle_DE))
dist_CE = np.sqrt((points['C'][0] - points['E'][0])**2 + (points['C'][1] - points['E'][1])**2)
return dist_CE - dist_CE_desired
def adjust_GJ(angle_GJ_deg, angle_BL_deg, dist_vars):
dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK = dist_vars
angle_BL = deg_to_rad(angle_BL_deg)
angle_GJ = deg_to_rad(angle_GJ_deg)
points = {
'A': (0, 0),
'B': (-60, -255),
'C': (-4, -669.5)
}
dist_BL = 2110
points['L'] = (points['B'][0] + dist_BL * np.cos(angle_BL), points['B'][1] + dist_BL * np.sin(angle_BL))
dist_BD = 1179.5
points['D'] = (points['B'][0] + dist_BD * np.cos(angle_BL), points['B'][1] + dist_BD * np.sin(angle_BL))
points['F'] = (points['B'][0] + dist_BF * np.cos(angle_BL), points['B'][1] + dist_BF * np.sin(angle_BL))
angle_DE = angle_BL + np.pi / 2
dist_DE = 30
points['E'] = (points['D'][0] + dist_DE * np.cos(angle_DE), points['D'][1] + dist_DE * np.sin(angle_DE))
if points['E'][1] <= points['D'][1]:
points['E'] = (points['D'][0] - dist_DE * np.cos(angle_DE), points['D'][1] - dist_DE * np.sin(angle_DE))
angle_GF = angle_BL + np.pi / 2
points['G'] = (points['F'][0] + dist_FG * np.cos(angle_GF), points['F'][1] + dist_FG * np.sin(angle_GF))
if points['G'][1] <= points['F'][1]:
points['G'] = (points['F'][0] - dist_FG * np.cos(angle_GF), points['F'][1] - dist_FG * np.sin(angle_GF))
points['J'] = (points['G'][0] + dist_GJ * np.cos(angle_GJ), points['G'][1] + dist_GJ * np.sin(angle_GJ))
if points['J'][1] <= points['G'][1]:
points['J'] = (points['G'][0] - dist_GJ * np.cos(angle_GJ), points['G'][1] - dist_GJ * np.sin(angle_GJ))
points['H'] = (points['G'][0] + dist_GH * np.cos(angle_GJ), points['G'][1] + dist_GH * np.sin(angle_GJ))
if points['H'][1] <= points['G'][1]:
points['H'] = (points['G'][0] - dist_GH * np.cos(angle_GJ), points['G'][1] - dist_GH * np.sin(angle_GJ))
angle_IH = angle_GJ + np.pi / 2
points['I'] = (points['H'][0] + dist_HI * np.cos(angle_IH), points['H'][1] + dist_HI * np.sin(angle_IH))
if points['I'][0] >= points['H'][0]:
points['I'] = (points['H'][0] - dist_HI * np.cos(angle_IH), points['H'][1] - dist_HI * np.sin(angle_IH))
def equations_K(vars, L, J):
x, y = vars
eq1 = (x - L[0])**2 + (y - L[1])**2 - 313.41**2
eq2 = (x - J[0])**2 + (y - J[1])**2 - dist_JK**2
return [eq1, eq2]
initial_guess1 = (points['L'][0] + 500, points['L'][1] + 500)
initial_guess2 = (points['L'][0] - 500, points['L'][1] - 500)
solution1 = fsolve(equations_K, initial_guess1, args=(points['L'], points['J']), xtol=1e-1)
solution2 = fsolve(equations_K, initial_guess2, args=(points['L'], points['J']), xtol=1e-1)
if np.allclose(solution1, solution2):
initial_guess1 = (points['L'][0] + 1000, points['L'][1] + 1000)
initial_guess2 = (points['L'][0] - 1000, points['L'][1] - 1000)
solution1 = fsolve(equations_K, initial_guess1, args=(points['L'], points['J']), xtol=1e-1)
solution2 = fsolve(equations_K, initial_guess2, args=(points['L'], points['J']), xtol=1e-1)
K = solution1 if solution1[1] > solution2[1] else solution2
points['K'] = K
dist_LM = 500
angle_ML = np.arctan2(points['L'][1] - K[1], points['L'][0] - K[0]) + deg_to_rad(109.638 + 180)
points['M'] = (points['L'][0] + dist_LM * np.cos(angle_ML), points['L'][1] + dist_LM * np.sin(angle_ML))
return abs(points['M'][1] - points['L'][1]), points
def compute_points(angle_BL_opt, angle_GJ_opt, dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK):
angle_BL = deg_to_rad(angle_BL_opt)
angle_GJ = deg_to_rad(angle_GJ_opt)
points = {
'A': (0, 0),
'B': (-60, -255),
'C': (-4, -669.5)
}
dist_BL = 2110
points['L'] = (points['B'][0] + dist_BL * np.cos(angle_BL), points['B'][1] + dist_BL * np.sin(angle_BL))
dist_BD = 1179.5
points['D'] = (points['B'][0] + dist_BD * np.cos(angle_BL), points['B'][1] + dist_BD * np.sin(angle_BL))
points['F'] = (points['B'][0] + dist_BF * np.cos(angle_BL), points['B'][1] + dist_BF * np.sin(angle_BL))
angle_DE = angle_BL + np.pi / 2
dist_DE = 30
points['E'] = (points['D'][0] + dist_DE * np.cos(angle_DE), points['D'][1] + dist_DE * np.sin(angle_DE))
if points['E'][1] <= points['D'][1]:
points['E'] = (points['D'][0] - dist_DE * np.cos(angle_DE), points['D'][1] - dist_DE * np.sin(angle_DE))
angle_GF = angle_BL + np.pi / 2
points['G'] = (points['F'][0] + dist_FG * np.cos(angle_GF), points['F'][1] + dist_FG * np.sin(angle_GF))
if points['G'][1] <= points['F'][1]:
points['G'] = (points['F'][0] - dist_FG * np.cos(angle_GF), points['F'][1] - dist_FG * np.sin(angle_GF))
points['J'] = (points['G'][0] + dist_GJ * np.cos(angle_GJ), points['G'][1] + dist_GJ * np.sin(angle_GJ))
if points['J'][1] <= points['G'][1]:
points['J'] = (points['G'][0] - dist_GJ * np.cos(angle_GJ), points['G'][1] - dist_GJ * np.sin(angle_GJ))
points['H'] = (points['G'][0] + dist_GH * np.cos(angle_GJ), points['G'][1] + dist_GH * np.sin(angle_GJ))
if points['H'][1] <= points['G'][1]:
points['H'] = (points['G'][0] - dist_GH * np.cos(angle_GJ), points['G'][1] - dist_GH * np.sin(angle_GJ))
angle_IH = angle_GJ + np.pi / 2
points['I'] = (points['H'][0] + dist_HI * np.cos(angle_IH), points['H'][1] + dist_HI * np.sin(angle_IH))
if points['I'][0] >= points['H'][0]:
points['I'] = (points['H'][0] - dist_HI * np.cos(angle_IH), points['H'][1] - dist_HI * np.sin(angle_IH))
def equations_K(vars, L, J):
x, y = vars
eq1 = (x - L[0])**2 + (y - L[1])**2 - 313.41**2
eq2 = (x - J[0])**2 + (y - J[1])**2 - dist_JK**2
return [eq1, eq2]
initial_guess1 = (points['L'][0] + 500, points['L'][1] + 500)
initial_guess2 = (points['L'][0] - 500, points['L'][1] - 500)
solution1 = fsolve(equations_K, initial_guess1, args=(points['L'], points['J']), xtol=1e-1)
solution2 = fsolve(equations_K, initial_guess2, args=(points['L'], points['J']), xtol=1e-1)
if np.allclose(solution1, solution2):
initial_guess1 = (points['L'][0] + 1000, points['L'][1] + 1000)
initial_guess2 = (points['L'][0] - 1000, points['L'][1] - 1000)
solution1 = fsolve(equations_K, initial_guess1, args=(points['L'], points['J']), xtol=1e-1)
solution2 = fsolve(equations_K, initial_guess2, args=(points['L'], points['J']), xtol=1e-1)
K = solution1 if solution1[1] > solution2[1] else solution2
points['K'] = K
dist_LM = 500
angle_ML = np.arctan2(points['L'][1] - K[1], points['L'][0] - K[0]) + deg_to_rad(109.638 + 180)
points['M'] = (points['L'][0] + dist_LM * np.cos(angle_ML), points['L'][1] + dist_LM * np.sin(angle_ML))
return points
def plot_points(points, title):
plt.figure(figsize=(10, 10))
for point, coord in points.items():
plt.plot(coord[0], coord[1], 'o', label=point)
plt.text(coord[0], coord[1], f" {point}", fontsize=12)
connections = [
('A', 'I'), ('A', 'B'),
('B', 'C'), ('B', 'L'),
('C', 'E'),
('D', 'E'),
('F', 'G'),
('G', 'J'),
('H', 'I'),
('J', 'K'),
('K', 'L'),
('L', 'M')
]
for start, end in connections:
if start in points and end in points:
plt.plot([points[start][0], points[end][0]], [points[start][1], points[end][1]], 'k-')
plt.xlabel('X')
plt.ylabel('Y')
plt.title(title)
plt.legend()
plt.grid(True)
plt.axis('equal')
plt.show()
def calculate_angle_between_vectors(p1, p2, p3, p4):
v1 = np.array([p2[0] - p1[0], p2[1] - p1[1]])
v2 = np.array([p4[0] - p3[0], p4[1] - p3[1]])
dot_product = np.dot(v1, v2)
magnitude_v1 = np.linalg.norm(v1)
magnitude_v2 = np.linalg.norm(v2)
cos_angle = dot_product / (magnitude_v1 * magnitude_v2)
angle = np.arccos(cos_angle)
return rad_to_deg(angle)
def minimize_diff_all(dist_vars, angle_BL_opt, angle_GJ_opt, dist_AI, initial_position=False):
dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK = dist_vars
points_initial = compute_points(angle_BL_opt, angle_GJ_opt, dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK)
diff_initial = abs(points_initial['M'][1] - points_initial['L'][1])
dist_CE_1300 = 1300
angle_BL_1300 = minimize_scalar(lambda x: np.abs(calculate_distance_CE(x, dist_CE_1300)), bounds=(0, 360), method='bounded').x
points_1300 = compute_points(angle_BL_1300, angle_GJ_opt, dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK)
diff_1300 = abs(points_1300['M'][1] - points_1300['L'][1])
dist_CE_1550 = 1550
angle_BL_1550 = minimize_scalar(lambda x: np.abs(calculate_distance_CE(x, dist_CE_1550)), bounds=(0, 360), method='bounded').x
points_1550 = compute_points(angle_BL_1550, angle_GJ_opt, dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK)
diff_1550 = abs(points_1550['M'][1] - points_1550['L'][1])
total_diff = diff_initial + diff_1300 + diff_1550
# Apply penalty for initial position if angle constraint is violated
if initial_position:
angle_JG_LA = calculate_angle_between_vectors(points_initial['J'], points_initial['G'], points_initial['L'], points_initial['A'])
if angle_JG_LA >= 160:
total_diff += 1e8 # Large penalty if the angle is not lower than 160 degrees
return total_diff
# Constraint function
def angle_constraint(dist_vars, angle_BL_opt, angle_GJ_opt):
dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK = dist_vars
points = compute_points(angle_BL_opt, angle_GJ_opt, dist_BF, dist_FG, dist_GJ, dist_GH, dist_HI, dist_JK)
angle_JG_LA = calculate_angle_between_vectors(points['J'], points['G'], points['L'], points['A'])
return 160 - angle_JG_LA # This should be greater than 0 for the constraint to be satisfied
dist_CE_desired = 1107
res_BL = minimize_scalar(lambda x: np.abs(calculate_distance_CE(x, dist_CE_desired)), bounds=(0, 360), method='bounded')
angle_BL_opt = res_BL.x
print(f"Optimal angle BL: {angle_BL_opt} degrees")
initial_guess = [1700, 0, 500, 300, 0, 700]
bounds = [(1500, 1900), (-50, 50), (400, 700), (100, 500), (-50, 50), (500, 900)]
# Optimize GJ with angle constraint applied for initial position
res_GJ = minimize_scalar(lambda x: adjust_GJ(x, angle_BL_opt, initial_guess)[0], bounds=(0, 360), method='bounded')
angle_GJ_opt = res_GJ.x
print(f"Optimal angle GJ: {angle_GJ_opt} degrees")
# Define the constraint as a NonlinearConstraint object
angle_constraint_obj = NonlinearConstraint(
lambda dist_vars: angle_constraint(dist_vars, angle_BL_opt, angle_GJ_opt),
0, # lower bound (constraint should be greater than 0)
np.inf # upper bound (no upper bound, we just want it to be greater than 0)
)
# Optimize distances with constraint applied only for the initial position
res = minimize(
lambda d: minimize_diff_all(d, angle_BL_opt, angle_GJ_opt, dist_AI, initial_position=True),
initial_guess,
bounds=bounds,
method='SLSQP',
constraints=[angle_constraint_obj]
)
optimized_distances = res.x
print(f"Optimized distances: {optimized_distances}")
points_initial = compute_points(angle_BL_opt, angle_GJ_opt, *optimized_distances)
angle_JG_LA_initial = calculate_angle_between_vectors(points_initial['J'], points_initial['G'], points_initial['L'], points_initial['A'])
print(f"Initial angle between JG and LA: {angle_JG_LA_initial} degrees")
if 'K' in points_initial:
print("Points after initial setup:")
for point, coord in points_initial.items():
print(f"{point}: {coord}")
print(f"Distances: BF={optimized_distances[0]}, FG={optimized_distances[1]}, GJ={optimized_distances[2]}, GH={optimized_distances[3]}, HI={optimized_distances[4]}, JK={optimized_distances[5]}")
plot_points(points_initial, f"Initial setup with CE = {dist_CE_desired} mm")
dist_CE_desired = 1300
angle_BL_1300 = minimize_scalar(lambda x: np.abs(calculate_distance_CE(x, dist_CE_desired)), bounds=(0, 360), method='bounded').x
points_1300 = compute_points(angle_BL_1300, angle_GJ_opt, *optimized_distances)
angle_JG_LA_1300 = calculate_angle_between_vectors(points_1300['J'], points_1300['G'], points_1300['L'], points_1300['A'])
print(f"Angle between JG and LA for CE=1300: {angle_JG_LA_1300} degrees")
if 'K' in points_1300:
print("Points after updating CE to 1300 mm:")
for point, coord in points_1300.items():
print(f"{point}: {coord}")
print(f"Distances: BF={optimized_distances[0]}, FG={optimized_distances[1]}, GJ={optimized_distances[2]}, GH={optimized_distances[3]}, HI={optimized_distances[4]}, JK={optimized_distances[5]}")
plot_points(points_1300, f"Updated setup with CE = 1300 mm and fixed AI = {dist_AI} mm")
dist_CE_desired = 1550
angle_BL_1550 = minimize_scalar(lambda x: np.abs(calculate_distance_CE(x, dist_CE_desired)), bounds=(0, 360), method='bounded').x
points_1550 = compute_points(angle_BL_1550, angle_GJ_opt, *optimized_distances)
angle_JG_LA_1550 = calculate_angle_between_vectors(points_1550['J'], points_1550['G'], points_1550['L'], points_1550['A'])
print(f"Angle between JG and LA for CE=1550: {angle_JG_LA_1550} degrees")
if 'K' in points_1550:
print("Points after updating CE to 1550 mm:")
for point, coord in points_1550.items():
print(f"{point}: {coord}")
print(f"Distances: BF={optimized_distances[0]}, FG={optimized_distances[1]}, GJ={optimized_distances[2]}, GH={optimized_distances[3]}, HI={optimized_distances[4]}, JK={optimized_distances[5]}")
plot_points(points_1550, f"Updated setup with CE = 1550 mm and fixed AI = {dist_AI} mm")
else:
print("Failed to compute valid points for initial setup.")
"""Part of the code that should allow to respect the angle constraint :
1.
def angle_constraint(dist_vars, angle_BL_opt, angle_GJ_opt):
angle_JG_LA = calculate_angle_between_vectors(points['J'], points['G'], points['L'], points['A'])
return 160 - angle_JG_LA # Should be greater than 0 to satisfy the constraint
2.
from scipy.optimize import NonlinearConstraint
angle_constraint_obj = NonlinearConstraint(
lambda dist_vars: angle_constraint(dist_vars, angle_BL_opt, angle_GJ_opt),
0, # Lower bound: angle should be less than 160 degrees (160 - angle >= 0)
np.inf # Upper bound: No upper bound, just need it greater than 0
)
3.
res = minimize(
lambda d: minimize_diff_all(d, angle_BL_opt, angle_GJ_opt, dist_AI, initial_position=True),
initial_guess,
bounds=bounds,
method='SLSQP',
constraints=[angle_constraint_obj] # Apply the angle constraint
)
"""
以下代码中，rad_to_deg角度转换不正确：
def deg_to_rad(angle_deg):
return angle_deg * np.pi / 180
def rad_to_deg(angle_rad):
return angle_rad * np.pi / 180
它正在做相反的转换：将度数转换为弧度。这个函数用在你的角度约束代码中，效果是它会一直认为角度小于0.05度。
更正了 rad_to_deg() 的代码：
def rad_to_deg(angle_rad):
return angle_rad / (np.pi / 180)