GeezerVison is showing github.com — Slycot CI check failing · Issue #1161 · python-control/python-control
Slycot CI check failing #1161DescriptionOver 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? Activitymurrayrm commented on Jun 25, 2025 murrayrm commented on Jun 26, 2025 ilayn commented on Jun 26, 2025 If I read the errors correctly, this is SLSQP output, we converted the F77 solver to C. murrayrm commented on Jun 28, 2025 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? ilayn commented on Jun 28, 2025 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? murrayrm commented on Jun 28, 2025 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. ilayn commented on Jun 30, 2025 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. ilayn commented on Jun 30, 2025 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. ilayn commented on Jun 30, 2025 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 costmurrayrm commented on Jun 30, 2025 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. murrayrm commented on Jul 1, 2025 (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...). ilayn commented on Jul 1, 2025 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. murrayrm commented on Jul 4, 2025 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... murrayrm commented on Jul 4, 2025 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). murrayrm commented on Jul 4, 2025 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. ilayn commented on Jul 5, 2025 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. murrayrm commented on Jul 5, 2025
MetadataMetadataDevelopmentNo branches or pull requestsIssue actions
|