Slycot CI check failing

Over the past few days, something changed and the Slycot CI check is now failing in GitHub actions. Here is an example:

https://github.com/python-control/python-control/actions/runs/15817454639/job/44617194878

The error is in control/tests/optimal_test and it looks like some of the optimization checks are no longer converging to the expected answer. Looks OK in the runs that use slycot from conda, so perhaps something different in the build dependencies?

NumPy 2.3 was released on 6 Jun 2025, so I'm guessing that is the culprit (somehow). Release notes:

https://numpy.org/doc/stable/release/2.3.0-notes.html

I've been able to recreate the problem locally. The problem is not NumPy 2.3, but SciPy 1.16.0:

https://docs.scipy.org/doc/scipy-1.16.0/release/1.16.0-notes.html

Working to figure out what the issue is ...

If I read the errors correctly, this is SLSQP output, we converted the F77 solver to C.

@ilayn

It looks like I'm getting a couple of different types of errors:

In some tests, I get res.success == False from scipy.optimize.minimize with message "Inequality constraints incompatible". I tried changing method to 'coblya' on one of these that I ran by hand and the error message went away.
In other tests, I get res.success == True but the optimal result doesn't match what it is expected to be. I tried setting tol=1e-12 and those error messages went away.

This gets rid of 13 of the 14 errors. Unfortunately, if I try to use coblya (or other solvers that can handle constraints) as the default method and I still get failures, mainly due to lack of convergence to the expected optimal values (not sure why yet).

@ilayn: Any insights into why the new C version might have differences in the way that constraints are handled and/or tolerances/convergence criteria? Is there a set of parameters I can pass via the options keyword that might recover the previous behavior?

I am a bit puzzled since in Ubuntu systems we have started seeing diverging results compared to the other architectures and I still don't have a good answer to these strange deviations. Is it also possible that the desired results are obtained by the F77 algorithm itself hence creating self-referencing results?

In any case, is it possible that we can isolate a simple case arguments that goes into scipy.optimize so I can try to emulate in a virtual box?

In the python-control context, here's an example of old versus new results:

import numpy as np
import scipy
import control as ct
import control.optimal as opt

# Define a (discrete-time) linear system with constraints
# Source: https://www.mpt3.org/UI/RegulationProblem

# LTI prediction model (discrete time)
sys = ct.ss([[1, 1], [0, 1]], [[1], [0.5]], np.eye(2), 0, 1)

# State and input constraints
constraints = [
    scipy.optimize.LinearConstraint(np.eye(3), [-5, -5, -1], [5, 5, 1]),
]

# Quadratic state and input penalty
Q = [[1, 0], [0, 1]]
R = [[1]]
cost = opt.quadratic_cost(sys, Q, R)

# Set up the optimal control problem
time = np.arange(0, 5, 1)
x0 = [4, 0]

# Retrieve the full open-loop predictions
res = opt.solve_optimal_trajectory(
    sys, time, x0, cost, constraints, squeeze=True,
    trajectory_method='shooting', 
    terminal_cost=cost)     # include to match MPT3 formulation

print("\nMPT inputs: [-1, -1, 0.1393, 0.3361, -5.204e-16]")
print("python-control: ", res.inputs)

print("\nMPT cost: 32.4898")
print("python-control: ", res.cost)

Results with SciPy 1.15.1:

Summary statistics:
* Cost function calls: 48
* Constraint calls: 63
* System simulations: 103
* Final cost: 32.48975410819071

MPT inputs: [-1, -1, 0.1393, 0.3361, -5.204e-16]
python-control:  [-1.         -1.          0.13929139  0.3361193   0.        ]

MPT cost: 32.4898
python-control:  32.48975410819071

The results from the MPT toolbox were created under MATLAB (though might still be using F77 SLSQP under the hood?).

Results with SciPy 1.16.0:

Summary statistics:
* Cost function calls: 56
* Constraint calls: 69
* System simulations: 89
* Final cost: 34.08095225022106

MPT inputs: [-1, -1, 0.1393, 0.3361, -5.204e-16]
python-control:  [-1.         -1.          0.21376303  0.78257832  1.        ]

MPT cost: 32.4898
python-control:  34.08095225022106

Since the cost is greater in the new case, it would seem that at the least it has different convergence criteria (or multiple minima, but I don't think that is the case here).

When I get a chance, I'll see if I can "unwrap" the optimization that python-control is setting up and come up with a simple case for the arguments that go directly into scipy.optimize.minimize.

Thank you for the extra information, will be very useful. I'll also try to build SciPy with logging so we can catch it on the other side too.

I just did a very crude printf() comparison. Even in the beginning when the function evaluation is requested at the point starts to disagree already

This is 1.15.1 output

 GS IS   -119.99999999097514     
 CONVERGENCE CHECK 170   235.99999999119444        0.0000000000000000        9.9999999999999995E-007
Mode 1: received x: [-1. -1. -1. -1.  0.]
Resulting function value: 47.99999999678176

and in the new version this gives

GS IS    -120.00000000086266
CONVERGENCE CHECK 170: h1 =     236.00000000084182, h2 =                      0, acc = 9.9999999999999995e-07
Mode 1: received x: [-1. -1. -1. -1.  0.]
Resulting function value: 84.00000000004164

Does this look meaningful to you? I could not dive into the problem setup so I can't quite see what the function is but it seems like the function evaluation is not the same between two versions. This might be completely SciPy's fault in its clamping behavior change but any feedback would be really nice.

From the inspect module I can see the function is the same to both versions so I'm not sure I can explain why they are giving different results yet

def _cost_function(self, coeffs):
    if self.log:
        start_time = time.process_time()
        logging.info("_cost_function called at: %g", start_time)

    # Compute the states and inputs
    states, inputs = self._compute_states_inputs(coeffs)

    # Trajectory cost
    if ct.isctime(self.system):
        # Evaluate the costs
        costs = [self.integral_cost(states[:, i], inputs[:, i]) for
                    i in range(self.timepts.size)]

        # Compute the time intervals
        dt = np.diff(self.timepts)

        # Integrate the cost
        costs = np.array(costs)
        # Approximate the integral using trapezoidal rule
        cost = np.sum(0.5 * (costs[:-1] + costs[1:]) * dt)

    else:
        # Sum the integral cost over the time (second) indices
        # cost += self.integral_cost(states[:,i], inputs[:,i])
        cost = sum(map(
            self.integral_cost, states[:, :-1].transpose(),
            inputs[:, :-1].transpose()))

    # Terminal cost
    if self.terminal_cost is not None:
        cost += self.terminal_cost(states[:, -1], inputs[:, -1])

    # Update statistics
    self.cost_evaluations += 1
    if self.log:
        stop_time = time.process_time()
        self.cost_process_time += stop_time - start_time
        logging.info(
            "_cost_function returning %g; elapsed time: %g",
            cost, stop_time - start_time)

    # Return the total cost for this input sequence
    return cost

I'll have a look. For this particular use case, what _cost_function expects as inputs is the values of input vector over time (SISO system => coeffs = (u[t], u[t+1], u[t+2], u[t+3], u[t+4]), as an ndarray). From those values, the system state over the horizon is generated via simulation (using odeint, I think).

The fact that the two versions are giving different costs is a great clue and it should help track down what is different.

(Mainly notes for me...)

I put in some statements in _cost_function() to show the inputs arguments (coeffs) as well as the states, inputs, costs. Here's what I get.

SciPy 1.15.1:

self.x=[4, 0]
coeffs=array([0., 0., 0., 0., 0.])
states=array([[4., 4., 4., 4., 4.],
       [0., 0., 0., 0., 0.]])
inputs=array([[0., 0., 0., 0., 0.]])
cost=80.0

## 5 iterations computing gradients (omitted)

self.x=[4, 0]
coeffs=array([-1., -1., -1., -1.,  0.])
states=array([[ 4. ,  3. ,  1.5, -0.5, -3. ],
       [ 0. , -0.5, -1. , -1.5, -2. ]])
inputs=array([[-1., -1., -1., -1.,  0.]])
cost=47.99999999854589

Scipy 1.16.0

self.x=[4, 0]
coeffs=array([0., 0., 0., 0., 0.])
states=array([[4., 4., 4., 4., 4.],
       [0., 0., 0., 0., 0.]])
inputs=array([[0., 0., 0., 0., 0.]])
cost=80.0

## 5 iterations computing gradients (omitted)

self.x=[4, 0]
coeffs=array([-1., -1., -1., -1.,  0.])
states=array([[4., 4., 4., 4., 4.],
       [0., 0., 0., 0., 0.]])
inputs=array([[-1., -1., -1., -1.,  0.]])
cost=83.9999999999217

So from this it looks like something is wrong with the computation of the states, inside the cost function. The initial state (first column) is the same in both cases and the inputs are the same, but for "SciPy 1.16.0" the input is somehow not getting applied.

Hard to imagine that is something in scipy.optimize.minimize, so it must be something else (maybe associated with SciPy 1.16.0 or maybe not...).

If SciPy is the only moving piece it is probably somewhere else but still in SciPy if I hazard a guess, odeint or solve_ivp.

Found the problem in optimal_test.py. Nothing to do with SciPy. Dumb coding error in saving the last value of a computation: needed to save a contents of an array, not just the reference to the array object. Not sure why it just started failing now...

Next problem comes up in lti_test.py, where ct.c2d() is throwing errors due to complex coefficients when creating a transfer function. Might be something different between NumPy 2.2.6 (old) and NumPy 2.3.1 (new).

Found the second set of errors. Looks like there is a difference in the way rounding was handled in NumPy 2.3.1 that wasn't showing up previously and arose when multiplying complex conjugates to get a real number.

PR coming next. I'll go ahead and merge once it comes through, since it is breaking CI and doesn't change functionality.

add bug fix associated with python-controlgh-1161

Ah OK, at least this time we dodged the bullet so that we don't need to play with SciPy version exceptions. Thank you for the extended investigation.

Fixed in #1163.

Navigation Menu

Repository navigation

Slycot CI check failing #1161

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Slycot CI check failing #1161

Description

Activity

murrayrm commented on Jun 25, 2025

murrayrm commented on Jun 26, 2025

ilayn commented on Jun 26, 2025

murrayrm commented on Jun 28, 2025

ilayn commented on Jun 28, 2025

murrayrm commented on Jun 28, 2025

ilayn commented on Jun 30, 2025

ilayn commented on Jun 30, 2025

ilayn commented on Jun 30, 2025

murrayrm commented on Jun 30, 2025

murrayrm commented on Jul 1, 2025

ilayn commented on Jul 1, 2025

murrayrm commented on Jul 4, 2025

murrayrm commented on Jul 4, 2025

murrayrm commented on Jul 4, 2025

ilayn commented on Jul 5, 2025

murrayrm commented on Jul 5, 2025

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions