Repository navigation
Slycot CI check failing #1161
Description
Activity
NumPy 2.3 was released on 6 Jun 2025, so I'm guessing that is the culprit (somehow). Release notes:
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.
It looks like I'm getting a couple of different types of errors:
-
In some tests, I get
res.success == Falsefromscipy.optimize.minimizewith message "Inequality constraints incompatible". I tried changingmethodto 'coblya' on one of these that I ran by hand and the error message went away. -
In other tests, I get
res.success == Truebut the optimal result doesn't match what it is expected to be. I tried settingtol=1e-12and 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
optionskeyword 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
desiredresults 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.optimizeso 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.48975410819071The 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.08095225022106Since 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 alreadyThis 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
inspectmodule 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 yetdef _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_functionexpects 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 (usingodeint, 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.99999999854589Scipy 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.9999999999217So 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,
odeintorsolve_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, wherect.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.
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.
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_testand it looks like some of the optimization checks are no longer converging to the expected answer. Looks OK in the runs that useslycotfrom conda, so perhaps something different in the build dependencies?