I have 4 datasets with the same length (701). The first represents the x values (station). The others the y values: one from a measurement (scannedCSLeft), and the rest are representing the permitted lowest and highest values for each station (loCSLeft and hiCSLeft). It looks like that: initial datasets
I wanted to generate with scipy.optimize.minimize a "best fit" polyline, which consists of just 10 points. The y values of the points are predefined by me, and the firts and last x values are also fixed. So the function has to define the rest 8 breakpoints x values.
I have 3 condition regarding the result: All points and all connecting lines must be between lowest and highest line (loCSLeft and hiCSLeft). Those lines, that connect breakpoints with different y values must have a slope between predefined constant values. And of course the result has to follow the measurement values as much as it can (the calculated error must be the lowest). The result should look like that (green line - but the slopes are not correct): result
I tried to solve this problem with scipy.optimize.minimize function, using 3 constraints, predefined bounds for each breakpoints, but I can't reach that point, when the result meets all 3 constraint.
I think the problem is with "constraint_slope_range" constraint, which should produce the correct slope for the non-horizontal sections. This constraint came from chatGPT, but I think the result should be a negative summarized number, to be harmonised with inequality constraint type (if I modify it to negative number the result become much worse, but I don't know why).
If anyone is interested, I can send you the datasets to be imported.
P.S. I'm just a beginner Python user, so there might be some messy parts in the script for you, but please get over it.
import os
import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as plt
########## IMPORTING PREPARED DATASETS ##########
file_path = os.path.dirname(os.path.abspath(__file__))
scannedCSLeft = np.genfromtxt(file_path+'/scannedCSLeft.csv')
loCSLeftMod = np.genfromtxt(file_path+'/loCSLeftMod.csv')
hiCSLeftMod = np.genfromtxt(file_path+'/hiCSLeftMod.csv')
# other inputs
station = np.linspace(0, 700, 701)
laneWidth = 3
min_slope = 30
max_slope = 50
# x coordinates of start and endpoint
xStart = station[0]
xEnd = station[-1]
# number of breakpoints, including the start and endpoint
noBreakPoints = 10
########## INITIAL BREAKPOINTS ##########
# initial x coordinates of breakpoints
x0 = np.linspace(xStart, xEnd, noBreakPoints)
# initial y coordinates of breakpoints - these are also the final values
y0 = (200, 200, -350, -350, 200, 200, 500, 500, 200, 200)
# array to be optimized: first and last item of x0 are fixed
opt = (x0[1:-1])
def get_xy_from_param_array(p):
# p contains the x coordinates of the not fixed breakpoints
x1 = p
x1 = np.concatenate(([xStart], x1, [xEnd]))
y1 = y0
return x1, y1
def multi_line_fit(p):
# querying the actual x-y coordinates of the breakpoints
x2, y2 = get_xy_from_param_array(p)
# interpolating the y coordinates in every round station
y3 = np.interp(station, x2, y2)
# the error value is the squared difference between the original measurment points (scannedCSLeft) and the interpolated points
error = np.sum((scannedCSLeft - y3)**2)
return error
def calculate_slope(p):
# querying the actual x-y coordinates of the breakpoints
x2, y2 = get_xy_from_param_array(p)
# calculating the slopes of each connecting line between the adjacent breakpoints
slopes = [abs(laneWidth * (y2[i+1] - y2[i]) / (x2[i+1] - x2[i])) for i in range(len(x2)-1) if laneWidth * (y2[i+1] - y2[i]) / (x2[i+1] - x2[i]) != 0]
return slopes
def constraint_min_y(p):
# every point and their connecting lines be above loCSLeft dataset
# querying the actual x-y coordinates of the breakpoints
x2, y2 = get_xy_from_param_array(p)
# interpolating the y coordinates in every round station
y3 = np.interp(station, x2, y2)
# difference between the interpolated points and the minimum points (loCSLeft)
diff = y3 - loCSLeftMod
# summarizing the negative differences, 'ineq' constraint will make it non-negative - every point and their connecting lines will be above loCSLeft dataset
sum_of_negatives = np.where(diff<0, diff, 0).sum()
return sum_of_negatives
def constraint_max_y(p):
# every point and their connecting lines be below hiCSLeft dataset
# querying the actual x-y coordinates of the breakpoints
x2, y2 = get_xy_from_param_array(p)
# interpolating the y coordinates in every round station
y3 = np.interp(station, x2, y2)
# difference between the maximum points (hiCSLeft) and the interpolated points
diff = hiCSLeftMod - y3
# summarizing the negative differences, 'ineq' constraint will make it non-negative - every point and their connecting lines will be below hiCSLeft dataset
sum_of_negatives = np.where(diff<0, diff, 0).sum()
return sum_of_negatives
def constraint_slope_range(p, min_slope, max_slope):
# Kényszerfüggvény, ami ellenőrzi, hogy az egyes egyenesek meredeksége a megadott határok között van-e
slopes = calculate_slope(p)
invalid_slopes = [slope for slope in slopes if slope < min_slope or slope > max_slope]
return np.sum(invalid_slopes)
constraints = [
{'type': 'ineq', 'fun': constraint_min_y},
{'type': 'ineq', 'fun': constraint_max_y},
{'type': 'ineq', 'fun': lambda p: constraint_slope_range(p, min_slope, max_slope)}
]
bounds = [(50, 150), (100, 200), (250, 350), (350, 450), (450, 550), (550, 600), (550, 650), (600, 700)]
res = minimize(fun=multi_line_fit, x0=opt, bounds=bounds, constraints=constraints)
y = y0
x = res.x
x = np.concatenate(([xStart], x, [xEnd]))
print(res.success)
print(res.message)
plt.plot(station, scannedCSLeft, linestyle='dotted')
plt.plot(station, loCSLeftMod, color='red', linestyle='dashed')
plt.plot(station, hiCSLeftMod, color='red', linestyle='dashed')
plt.plot(x, y, color='green', linewidth=3)
plt.show()