Add attribute to overwrite the velocity constraints

[tjstienstra]

References to other Issues or PRs

Brief description of what is fixed or changed

This PR adds the velocity_constraints attribute to System to allow users to overwrite the velocity constraints, which can be computed automatically in most cases. An example where one is required to overwrite the velocity constraints is when computing the noncontributing forces in the closing joint of a four-bar linkage.

Other comments

The velocity constraints are only stored directly when specifically set by the user. Otherwise, they are computed based on the holonomic and nonholonomic constraints. This is to prevent possible inconsistencies in the system as the same constraint can/should be stored as both a nonholonomic constraint as well as a velocity constraint.

Release Notes

• physics.mechanics
  • Implemented System.velocity_constraints to allow overwriting the automatically computed velocity constraints.

[sympy-bot]

✅

Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• physics.mechanics
  • Implemented System.velocity_constraints to allow overwriting the automatically computed velocity constraints. (#26374 by @tjstienstra)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed
This PR adds the `velocity_constraints` attribute to `System` to allow users to overwrite the velocity constraints, which can be computed automatically in most cases. An example where one is required to overwrite the velocity constraints is when computing the noncontributing forces in the closing joint of a four-bar linkage.

#### Other comments
The velocity constraints are only stored directly when specifically set by the user. Otherwise, they are computed based on the holonomic and nonholonomic constraints. This is to prevent possible inconsistencies in the system as the same constraint can/should be stored as both a nonholonomic constraint as well as a velocity constraint.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* physics.mechanics
  * Implemented System.velocity_constraints to allow overwriting the automatically computed velocity constraints.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[moorepants]

Does that value look right? Magnitude an direction? It should be the proportion of weight on one end of how ever it is configured?

[tjstienstra]

The four-bar linkage is oriented as a perfect rectangle with link2 and link4 oriented upwards. Therefore, half of the mass of link3, 0.5 * 2 * 5 = 5, and the full mass of link4, 1 * 5 = 5 rests on this joint. So a force of 9.81 * (5 + 5) = 98.1 acts downwards on the closing joint.

[moorepants]

This looks good to me and nice addition to the example. I'm not sure we have any example that shows how to do this for a closed loop any where. I left a couple small comments and maybe you can do a back of than calculation to see if that value is reasonable (if you haven't already).

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

| Change   | Before [6d28cd91]    | After [15db534f]    |   Ratio | Benchmark (Parameter)                                                |
|----------|----------------------|---------------------|---------|----------------------------------------------------------------------|
| -        | 70.1±0.8ms           | 44.7±0.5ms          |    0.64 | integrate.TimeIntegrationRisch02.time_doit(10)                       |
| -        | 68.0±0.7ms           | 43.2±0.3ms          |    0.64 | integrate.TimeIntegrationRisch02.time_doit_risch(10)                 |
| +        | 18.3±0.1μs           | 29.6±0.1μs          |    1.62 | integrate.TimeIntegrationRisch03.time_doit(1)                        |
| -        | 5.25±0.03ms          | 2.90±0.04ms         |    0.55 | logic.LogicSuite.time_load_file                                      |
| -        | 72.3±0.8ms           | 28.9±0.2ms          |    0.4  | polys.TimeGCD_GaussInt.time_op(1, 'dense')                           |
| -        | 26.0±0.3ms           | 17.2±0.09ms         |    0.66 | polys.TimeGCD_GaussInt.time_op(1, 'expr')                            |
| -        | 73.1±0.4ms           | 29.2±0.06ms         |    0.4  | polys.TimeGCD_GaussInt.time_op(1, 'sparse')                          |
| -        | 253±3ms              | 125±0.2ms           |    0.49 | polys.TimeGCD_GaussInt.time_op(2, 'dense')                           |
| -        | 254±2ms              | 125±0.8ms           |    0.49 | polys.TimeGCD_GaussInt.time_op(2, 'sparse')                          |
| -        | 659±3ms              | 372±0.9ms           |    0.57 | polys.TimeGCD_GaussInt.time_op(3, 'dense')                           |
| -        | 652±4ms              | 377±1ms             |    0.58 | polys.TimeGCD_GaussInt.time_op(3, 'sparse')                          |
| -        | 497±3μs              | 291±2μs             |    0.59 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(1, 'dense')            |
| -        | 1.78±0.01ms          | 1.06±0.01ms         |    0.6  | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(2, 'dense')            |
| -        | 5.80±0.04ms          | 3.16±0.02ms         |    0.54 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(3, 'dense')            |
| -        | 446±3μs              | 231±0.8μs           |    0.52 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(1, 'dense')               |
| -        | 1.48±0.01ms          | 677±5μs             |    0.46 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(2, 'dense')               |
| -        | 4.90±0.03ms          | 1.65±0.01ms         |    0.34 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(3, 'dense')               |
| -        | 376±0.5μs            | 206±1μs             |    0.55 | polys.TimeGCD_SparseGCDHighDegree.time_op(1, 'dense')                |
| -        | 2.42±0.01ms          | 1.24±0.01ms         |    0.51 | polys.TimeGCD_SparseGCDHighDegree.time_op(3, 'dense')                |
| -        | 10.1±0.02ms          | 4.47±0.01ms         |    0.44 | polys.TimeGCD_SparseGCDHighDegree.time_op(5, 'dense')                |
| -        | 354±5μs              | 171±0.7μs           |    0.48 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(1, 'dense')            |
| -        | 2.47±0.01ms          | 902±8μs             |    0.37 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(3, 'dense')            |
| -        | 9.51±0.02ms          | 2.64±0.02ms         |    0.28 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(5, 'dense')            |
| -        | 1.02±0.01ms          | 430±3μs             |    0.42 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'dense')           |
| -        | 1.73±0.02ms          | 504±3μs             |    0.29 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')          |
| -        | 5.87±0.06ms          | 1.80±0.02ms         |    0.31 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'dense')           |
| -        | 8.43±0.09ms          | 1.51±0.01ms         |    0.18 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')          |
| -        | 283±1μs              | 65.0±0.3μs          |    0.23 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')             |
| -        | 3.44±0.04ms          | 393±1μs             |    0.11 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'dense')              |
| -        | 4.00±0.02ms          | 279±0.8μs           |    0.07 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')             |
| -        | 7.08±0.03ms          | 1.28±0ms            |    0.18 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'dense')              |
| -        | 8.77±0.04ms          | 846±5μs             |    0.1  | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')             |
| -        | 4.99±0ms             | 3.01±0.03ms         |    0.6  | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| -        | 11.9±0.1ms           | 6.66±0.03ms         |    0.56 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'dense')  |
| -        | 22.3±0.08ms          | 9.04±0.03ms         |    0.4  | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| -        | 5.22±0.02ms          | 870±4μs             |    0.17 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')    |
| -        | 12.4±0.04ms          | 7.13±0.06ms         |    0.57 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')    |
| -        | 101±0.7ms            | 25.8±0.08ms         |    0.25 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'dense')     |
| -        | 166±0.8ms            | 54.2±0.3ms          |    0.33 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')    |
| -        | 177±3μs              | 112±0.5μs           |    0.64 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'dense')      |
| -        | 355±2μs              | 217±2μs             |    0.61 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'sparse')     |
| -        | 4.24±0.03ms          | 833±2μs             |    0.2  | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'dense')      |
| -        | 5.28±0.04ms          | 386±2μs             |    0.07 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')     |
| -        | 19.9±0.3ms           | 2.80±0.03ms         |    0.14 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'dense')      |
| -        | 22.7±0.1ms           | 626±0.8μs           |    0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')     |
| -        | 482±5μs              | 138±0.8μs           |    0.29 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| -        | 4.68±0.04ms          | 628±3μs             |    0.13 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'dense')  |
| -        | 5.35±0.02ms          | 140±0.4μs           |    0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| -        | 13.4±0.3ms           | 1.31±0.02ms         |    0.1  | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'dense')  |
| -        | 13.9±0.06ms          | 142±2μs             |    0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| -        | 132±1μs              | 75.0±1μs            |    0.57 | solve.TimeMatrixOperations.time_rref(3, 0)                           |
| -        | 247±2μs              | 87.7±0.5μs          |    0.35 | solve.TimeMatrixOperations.time_rref(4, 0)                           |
| -        | 24.5±0.08ms          | 10.3±0.05ms         |    0.42 | solve.TimeSolveLinSys189x49.time_solve_lin_sys                       |
| -        | 28.6±0.3ms           | 15.7±0.1ms          |    0.55 | solve.TimeSparseSystem.time_linsolve_Aaug(20)                        |
| -        | 55.1±0.5ms           | 25.3±0.09ms         |    0.46 | solve.TimeSparseSystem.time_linsolve_Aaug(30)                        |
| -        | 28.3±0.2ms           | 15.5±0.2ms          |    0.55 | solve.TimeSparseSystem.time_linsolve_Ab(20)                          |
| -        | 54.8±0.2ms           | 25.3±0.3ms          |    0.46 | solve.TimeSparseSystem.time_linsolve_Ab(30)                          |

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).