Skip to content

Friction Oscillator

System Description

Self-excited friction oscillator with Stribeck friction characteristic:

\[ \begin{aligned} \dot{x} &= v \\ \dot{v} &= -2\xi v - \mu(v_{\text{rel}}) \cdot \text{sgn}(v_{\text{rel}}) \end{aligned} \]

where \(v_{\text{rel}} = v - v_d\) and the friction characteristic is:

\[ \mu(v_{\text{rel}}) = \mu_d + (\mu_{sd} - 1)\mu_d \exp\left(-\frac{|v_{\text{rel}}|}{v_0}\right) + \mu_v |v_{\text{rel}}| \]

System Parameters

Parameter Symbol Value
Driving velocity \(v_d\) 1.5
Damping ratio \(\xi\) 0.05
Static/dynamic ratio \(\mu_{sd}\) 2.0
Dynamic friction \(\mu_d\) 0.5
Viscous friction \(\mu_v\) 0.0
Reference velocity \(v_0\) 0.5

Sampling

  • Dimension: \(D = 2\)
  • Sample size: \(N = 5000\)
  • Distribution: \(\rho\) = Uniform
  • Region of interest: \(\mathcal{Q}(x, v) : [-2, 2] \times [0, 2]\)

Solver

Setting Value
Method Dopri5 (Diffrax)
Time span \([0, 500]\)
Steps 500 (\(f_s\) = 1 Hz)
Relative tolerance 1e-08
Absolute tolerance 1e-06

Implementation detail

The JAX implementation uses a \(\tanh\) approximation for the sign function with \(k_{\text{smooth}} = 200\), while MATLAB uses hard switching.

Feature Extraction

Maximum absolute velocity after transient:

  • States: \(v\) (state 1)
  • Formula: \(\max(|v|)\), threshold = 0.2 (FP if \(\leq 0.2\), LC otherwise)
  • Transient cutoff: \(t^* = 400.0\)

Clustering

  • Method: k-NN (k=1)
  • Template ICs:
  • FP: \([1.0, 1.0]\) — Fixed point (stable fixed point of steady sliding)
  • LC: \([2.0, 2.0]\) — Limit cycle (stick-slip oscillation)

Reproduction Code

Setup

def setup_friction_system() -> SetupProperties:
    n = 5000

    device = "cuda" if torch.cuda.is_available() else "cpu"
    print(f"Setting up friction system on device: {device}")

    params: FrictionParams = {
        "v_d": 1.5,  # Driving velocity
        "xi": 0.05,  # Damping ratio
        "musd": 2.0,  # Ratio static to dynamic friction coefficient
        "mud": 0.5,  # Dynamic coefficient of friction
        "muv": 0.0,  # Linear strengthening parameter
        "v0": 0.5,  # Reference velocity for exponential decay
    }

    ode_system = FrictionJaxODE(params)

    sampler = UniformRandomSampler(
        min_limits=[-2.0, 0.0],
        max_limits=[2.0, 2.0],
        device=device,
    )

    solver = JaxSolver(
        t_span=(0, 500),
        t_steps=500,
        t_eval=(400.0, 500.0),
        device=device,
        rtol=1e-8,
        atol=1e-6,
        cache_dir=".pybasin_cache/friction",
    )

    feature_extractor = FrictionFeatureExtractor(time_steady=400)

    classifier_initial_conditions = [
        [1.0, 1.0],
        [2.0, 2.0],
    ]

    classifier_labels = ["FP", "LC"]

    knn = KNeighborsClassifier(n_neighbors=1)

    template_integrator = TemplateIntegrator(
        template_y0=classifier_initial_conditions,
        labels=classifier_labels,
        ode_params=params,
    )

    return {
        "n": n,
        "ode_system": ode_system,
        "sampler": sampler,
        "solver": solver,
        "feature_extractor": feature_extractor,
        "estimator": knn,
        "template_integrator": template_integrator,
    }

Main Estimation

def main():
    props = setup_friction_system()

    bse = BasinStabilityEstimator(
        n=props["n"],
        ode_system=props["ode_system"],
        sampler=props["sampler"],
        solver=props.get("solver"),
        feature_extractor=props.get("feature_extractor"),
        predictor=props.get("estimator"),
        template_integrator=props.get("template_integrator"),
        output_dir="results_friction",
        feature_selector=None,
    )

    bse.run()

    return bse

Case 1: Baseline Results

Comparison with MATLAB bSTAB

Overall Classification Quality:

  • Matthews Correlation Coefficient: 1.0000
Attractor pybasin BS ± SE bSTAB BS ± SE
FP 0.31200 ± 0.00655 0.31200 ± 0.00655
LC 0.68800 ± 0.00655 0.68800 ± 0.00655

Visualizations

Basin Stability

Basin Stability

State Space

State Space

Feature Space

Feature Space

Template Trajectories

Trajectories

Case 2: v_d Parameter Sweep

Comparison with MATLAB bSTAB

Average MCC = 0.9991

The average excludes cases where there is only a single attractor and the basin stability values are the same since MCC is 0 for single class cases, and would therefore drop the average.

Parameter Attractor pybasin BS ± SE bSTAB BS ± SE MCC
0.8 LC 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
0.875 LC 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
0.95 LC 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
1.025 LC 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
1.1 LC 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
1.175 FP 0.01280 ± 0.00159 0.01260 ± 0.00158 0.9921
LC 0.98720 ± 0.00159 0.98740 ± 0.00158
1.25 FP 0.08380 ± 0.00392 0.08380 ± 0.00392 1.0000
LC 0.91620 ± 0.00392 0.91620 ± 0.00392
1.325 FP 0.13780 ± 0.00487 0.13780 ± 0.00487 1.0000
LC 0.86220 ± 0.00487 0.86220 ± 0.00487
1.4 FP 0.21100 ± 0.00577 0.21100 ± 0.00577 1.0000
LC 0.78900 ± 0.00577 0.78900 ± 0.00577
1.475 FP 0.28240 ± 0.00637 0.28240 ± 0.00637 1.0000
LC 0.71760 ± 0.00637 0.71760 ± 0.00637
1.55 FP 0.36360 ± 0.00680 0.36360 ± 0.00680 1.0000
LC 0.63640 ± 0.00680 0.63640 ± 0.00680
1.625 FP 0.43360 ± 0.00701 0.43360 ± 0.00701 1.0000
LC 0.56640 ± 0.00701 0.56640 ± 0.00701
1.7 FP 0.49120 ± 0.00707 0.49120 ± 0.00707 1.0000
LC 0.50880 ± 0.00707 0.50880 ± 0.00707
1.775 FP 0.55520 ± 0.00703 0.55520 ± 0.00703 1.0000
LC 0.44480 ± 0.00703 0.44480 ± 0.00703
1.85 FP 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
1.925 FP 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
2 FP 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
2.075 FP 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
2.15 FP 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000
2.225 FP 1.00000 ± 0.00000 1.00000 ± 0.00000 0.0000

Visualizations

Basin Stability Variation

Basin Stability Variation

Bifurcation Diagram

Bifurcation Diagram