Resolving trade off between smoothing via polynomial fitting and simplify-like algorithms (do not eliminate vertices)

I am struggling for a few weeks to use Python to smooth various linestrings and bumped into a conundrum. There seems to be a trade-off between using techniques that apply parametric curve fitting and those that apply Non-parametric algorithms like shapely simplify. I have attempted to use the following parametric strategies:

I used numpy polyfit and polyval

coefs_xy = np.polyfit(xcoords, ycoords, 4)
fitted_XY = np.polyval(coefs_xy, xcoords)

that fit polynomial functions to x and y coordinates. And used scipy UnivariateSpline that does a similar process.

spl = UnivariateSpline(xaxis, yaxis, s=.01, k=5)

They work fine if linestrings are well-behaved, i.e., if an x coordinate has only one corresponding y coordinate. This is violated when a linestring is a vertical arch, for example. In this case, the algorithm interprets two y coordinates as a statistical error around a mean and tries to find a ‘mean’ y coordinate between the two points. Or, in the scipy case, it starts zigzagging between y coordinates. Of course, the smoothing goes to the drain. I have even attempted to invert x with y coordinates and pick the one that has the smallest deviation wrt original linestring. However, there is always an exception that causes problems.

Alternatively, I have applied shapely simplify. It does solve the problem above. It may smooth corners locally, but looking to the whole line, I can still see vertices, especially when you have an almost straight line then a small angle. This is the code I used:

def interpolateBySimplify(lineCoords, crs):
    lineToBeSimplified = LineString(lineCoords)
    denseLine = np.linspace(0, lineToBeSimplified.length, 1000)
    tempVertList = []
    for vert in denseLine:
        tempPoint = lineToBeSimplified.interpolate(vert)
        tempVertList.append(tempPoint)
    interpolatedLine = []
    for i in range(len(tempVertList)):
        interpolatedLine.append(tempVertList[i].coords[0])
    lineToBeSimplified = LineString(interpolatedLine)
    gdf = gpd.GeoDataFrame(columns =['id', 'geometry'], crs=crs)
    new_row = {'id':0, 'geometry':lineToBeSimplified}
    gdf = gdf.append(new_row, ignore_index=True)
    tolerance = 3
    simplified = gdf.simplify(tolerance, preserve_topology=False)
    simplified.buffer(200, join_style=1).buffer(-200, join_style=1)
    finalLine_interpolated = {'id': [1], 'geometry':[LineString(lineToBeSimplified)]}
    return finalLine_interpolated, interpolatedLine

Even with a very large buffer (this is UTM coordinates), I still can see vertices like those in the figure below

Have you come across this problem?

It may be the case I am missing some detail in my approach.