MohrCoulomb-SimplexFit

A reproducible Nelder–Mead framework for inverse identification of Mohr–Coulomb shear-strength parameters from triaxial failure data

MohrCoulomb-SimplexFit is a research-oriented command-line application and Python library for estimating cohesion, $c$, and friction angle, $\varphi$, from conventional triaxial compression failure stresses. The program maps each observed principal-stress pair to a Mohr circle, formulates parameter identification as a constrained geometric tangency problem, and solves the resulting nonlinear least-squares objective with a native Nelder–Mead simplex implementation.

The software is designed around four principles: a mathematically explicit objective, physically admissible parameters, deterministic computation, and auditable outputs. SciPy is used only as a development-time reference and is not required at runtime.

Scientific scope

Given $n$ triaxial failure observations

$$ \mathcal{D}={(\sigma_{3i},\sigma_{1i})}_{i=1}^{n}, $$

the program identifies the linear Mohr–Coulomb envelope

$$ \tau=c+\sigma\tan\varphi. $$

Version 0.1 adopts the following assumptions:

• compression is positive and $0\leq\sigma_{3i}\leq\sigma_{1i}$;
• every input value uses one consistent stress unit;
• at least three observations and two distinct stress states are available;
• the failure envelope is linear over the investigated stress range;
• all observations receive equal weight;
• the admissible parameter domain is $c\geq0$ and $0\leq\varphi<90^\circ$.

Here, $\sigma_1$ is the major principal stress, $\sigma_3$ is the minor principal stress, $\sigma$ is normal stress, and $\tau$ is shear stress. Total or effective stresses may be used, but they must not be mixed within one analysis.

Mathematical formulation and derivation

1. From principal stresses to a Mohr circle

For observation $i$, the two-dimensional stress-transformation equations are

$$ \sigma_n(\theta)=\frac{\sigma_{1i}+\sigma_{3i}}{2} +\frac{\sigma_{1i}-\sigma_{3i}}{2}\cos 2\theta, $$

$$ \tau_n(\theta)=\frac{\sigma_{1i}-\sigma_{3i}}{2}\sin 2\theta. $$

Defining

$$ O_i=\frac{\sigma_{1i}+\sigma_{3i}}{2}, \qquad R_i=\frac{\sigma_{1i}-\sigma_{3i}}{2}, $$

and eliminating $\theta$ gives

$$ (\sigma_n-O_i)^2+\tau_n^2=R_i^2. $$

Thus each triaxial test is represented by a circle centered at $(O_i,0)$ with radius $R_i$.

2. Geometric tangency condition

The strength envelope can be written in implicit line form as

$$ \sigma\tan\varphi-\tau+c=0. $$

The perpendicular distance from the center $(O_i,0)$ to this line is

$$ D_i(c,\varphi) =\frac{|O_i\tan\varphi+c|}{\sqrt{\tan^2\varphi+1}}. $$

For the adopted physical domain, $O_i\geq0$, $c\geq0$, and $0\leq\varphi<90^\circ$. The absolute value can therefore be removed. Using $\sqrt{1+\tan^2\varphi}=1/\cos\varphi$,

$$ D_i(c,\varphi)=O_i\sin\varphi+c\cos\varphi. $$

Exact tangency requires the center-to-line distance to equal the circle radius:

$$ R_i=O_i\sin\varphi+c\cos\varphi. $$

The signed geometric tangency residual used by the software is consequently

$$ e_i(c,\varphi) =R_i-O_i\sin\varphi-c\cos\varphi. $$

A positive residual means that the candidate line intersects the ideal tangency distance of that circle; a negative residual means that the line lies farther from the center than the radius.

3. Connection to the conventional triaxial failure equation

Substituting the definitions of $O_i$ and $R_i$ into the exact tangency condition gives

$$ \frac{\sigma_{1i}-\sigma_{3i}}{2} =\frac{\sigma_{1i}+\sigma_{3i}}{2}\sin\varphi+c\cos\varphi. $$

After collecting the principal-stress terms,

$$ \sigma_{1i}(1-\sin\varphi) =\sigma_{3i}(1+\sin\varphi)+2c\cos\varphi, $$

and hence

$$ \sigma_{1i} =\sigma_{3i}\frac{1+\sin\varphi}{1-\sin\varphi} +\frac{2c\cos\varphi}{1-\sin\varphi}. $$

Using the half-angle identities yields the familiar form

$$ \sigma_{1i} =\sigma_{3i}\tan^2(45^\circ+\frac{\varphi}{2}) +2c\tan(45^\circ+\frac{\varphi}{2}). $$

This establishes the equivalence between Mohr-circle tangency and the standard triaxial Mohr–Coulomb failure relationship. The implementation minimizes geometric tangency residuals in the $(\sigma,\tau)$ plane rather than vertical residuals in $\sigma_1$.

4. Inverse problem and error measures

For imperfect experimental data, one line will not generally be tangent to every circle. The unweighted inverse problem is

$$ (\hat c,\hat\varphi) =\arg\min_{c,\varphi} J(c,\varphi), \qquad J(c,\varphi)=\sum_{i=1}^{n}e_i(c,\varphi)^2, $$

subject to

$$ c\geq0, \qquad 0\leq\varphi<\frac{\pi}{2}. $$

The reported objective has squared-stress units. A directly interpretable summary is

$$ \mathrm{RMSE} =\sqrt{\frac{J(\hat c,\hat\varphi)}{n}}, $$

which has the same unit as the input stresses. Residual signs and magnitudes remain available for every observation in residuals.csv; RMSE should be interpreted together with the stress scale and the residual pattern, not as a universal goodness-of-fit threshold.

5. Initialization and numerical scaling

The tangency equation is linear after introducing

$$ a=\sin\varphi, \qquad b=c\cos\varphi, \qquad R_i=aO_i+b. $$

The program first solves the auxiliary least-squares system

$$ \begin{bmatrix} O_1 & 1\\ \vdots & \vdots\\ O_n & 1 \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} \approx \begin{bmatrix} R_1 \\ \vdots \\ R_n \end{bmatrix}. $$

The estimate is projected onto $0\leq a&lt;1$, then transformed through

$$ \varphi_0=\arcsin a, \qquad c_0=\max(0,\frac{b}{\cos\varphi_0}). $$

This deterministic initializer reduces sensitivity to a poor arbitrary starting simplex. To limit scale imbalance between cohesion and angle, the optimization coordinates are

$$ \mathbf{x}=[\frac{c}{S},\ \varphi]^T, \qquad S=\max_i\sigma_{1i}, $$

and the minimized internal objective is $\widetilde J=J/S^2$. Reported parameters, residuals, $J$, and RMSE are converted back to the original input unit.

6. Native Nelder–Mead simplex search

For the two-parameter vector $\mathbf{x}$, the simplex has three vertices. At each iteration the vertices are ordered so that

$$ f(\mathbf{x}_1)\leq f(\mathbf{x}_2)\leq f(\mathbf{x}_3), $$

and the centroid excluding the worst vertex is

$$ \bar{\mathbf{x}}=\frac{\mathbf{x}_1+\mathbf{x}_2}{2}. $$

The implementation applies the standard operations:

Operation	Trial point	Coefficient
Reflection	$\mathbf{x}_r=\bar{\mathbf{x}}+\alpha(\bar{\mathbf{x}}-\mathbf{x}_3)$	$\alpha=1$
Expansion	$\mathbf{x}_e=\bar{\mathbf{x}}+\gamma(\mathbf{x}_r-\bar{\mathbf{x}})$	$\gamma=2$
Outside contraction	$\mathbf{x}_{oc}=\bar{\mathbf{x}}+\rho(\mathbf{x}_r-\bar{\mathbf{x}})$	$\rho=0.5$
Inside contraction	$\mathbf{x}_{ic}=\bar{\mathbf{x}}-\rho(\bar{\mathbf{x}}-\mathbf{x}_3)$	$\rho=0.5$
Shrink	$\mathbf{x}_i\leftarrow\mathbf{x}_1+\eta(\mathbf{x}_i-\mathbf{x}_1)$	$\eta=0.5$

Every trial vertex is projected onto the admissible domain. Convergence requires both the maximum simplex coordinate span and maximum objective-value span to fall below their configured tolerances. The default iteration limit is 2000; the result records the iteration count, function-evaluation count, convergence state, and termination reason.

Reproducible example

The repository includes three triaxial failure observations in examples/triaxial_example.csv:

Test	$\sigma_3$ (kPa)	$\sigma_1$ (kPa)
test_1	100.1	594.9357
test_2	200.3	1147.3320
test_3	301.2	1574.3210

The deterministic v0.1 analysis gives

$$ \hat c=28.115565\ \mathrm{kPa}, \qquad \hat\varphi=41.345419^\circ, $$

with $J=79.961597\ \mathrm{kPa}^2$ and $\mathrm{RMSE}=5.162738\ \mathrm{kPa}$. The figure below is generated directly from the tracked CSV and the public fitting API; it combines the fitted envelope, local objective landscape, and observation-level diagnostics.

Regenerate both scientific figures with

python scripts/generate_readme_figures.py

Installation

Python 3.10 or newer is required. For a conventional installation:

python -m venv .venv
python -m pip install .

Install development tools, including the SciPy reference implementation, with

python -m pip install -e ".[dev]"

Routine analysis workflow

Interactive launcher

Place custom CSV files in data/input/, then run the launcher for your platform.

Windows PowerShell:

.\run_fit.ps1

Linux, macOS, or Git Bash:

chmod +x run_fit.sh
./run_fit.sh

The launcher creates .venv when required, installs the project, lists available CSV files, asks for a stress unit, and writes a timestamped analysis under results/. Non-interactive execution is also available:

.\run_fit.ps1 -InputFile data\input\triaxial_example.csv -StressUnit kPa -OutputName my-run

./run_fit.sh --input data/input/triaxial_example.csv --unit kPa --output-name my-run

Direct CLI

mohrcoulomb-simplexfit fit examples/triaxial_example.csv \
  --stress-unit kPa \
  --output-dir results/example

Small datasets may be entered without creating a file:

mohrcoulomb-simplexfit fit \
  --pair 100.1 594.9357 \
  --pair 200.3 1147.332 \
  --pair 301.2 1574.321 \
  --stress-unit kPa \
  --output-dir results/inline-example

Use mohrcoulomb-simplexfit fit --help for initial-value and tolerance controls. Non-empty output directories are protected; --overwrite replaces only the three standard result files and preserves unrelated files.

Input and output contract

CSV input requires sigma3 and sigma1; label is optional and additional columns are ignored.

label,sigma3,sigma1
specimen_01,100,500
specimen_02,200,800
specimen_03,300,1100

The --stress-unit value is metadata, not a conversion instruction. Each successful analysis creates:

Artifact	Purpose
result.json	Versioned machine-readable parameters, input metadata, and optimizer diagnostics
residuals.csv	Observation-level circles, distances, signed residuals, and squared residuals
mohr_coulomb_fit.png	Equal-axis Mohr circles and fitted upper failure envelope

Validation errors return process exit code 2; numerical convergence or output failures return 1; successful analyses return 0.

Python API

from mohrcoulomb_simplexfit import FitOptions, StressPair, fit

observations = [
    StressPair(100.1, 594.9357, "test_1"),
    StressPair(200.3, 1147.332, "test_2"),
    StressPair(301.2, 1574.321, "test_3"),
]

result = fit(observations, options=FitOptions(max_iterations=2000))
print(result.cohesion)
print(result.friction_angle_degrees)
print(result.rmse)

Invalid data or controls raise ValidationError; failure to meet numerical convergence criteria raises ConvergenceError.

Verification and software quality

The automated suite covers analytical circle geometry, residual evaluation, synthetic exact parameter recovery, stress scaling, noisy data, physical boundaries, optimizer failure modes, CSV validation, CLI exit codes, and output artifacts. The native result is cross-validated against scipy.optimize.minimize(method="Nelder-Mead") on deterministic datasets.

ruff check .
ruff format --check .
pytest
python -m build

The present suite contains 24 tests and runs on Python 3.10 and 3.12 in CI.

Limitations and responsible interpretation

• The model is deterministic and does not produce confidence or prediction intervals.
• The objective is unweighted; heteroscedastic measurement uncertainty is not represented.
• Linear Mohr–Coulomb behavior is assumed over the supplied confining-stress range.
• Pore-pressure corrections and effective-stress conversion must be completed before fitting.
• Correlated replicates, censored tests, tensile states, nonlinear envelopes, and model-selection diagnostics are outside the v0.1 scope.
• A converged numerical optimum is not, by itself, evidence that the constitutive model is adequate. Inspect the residual table, fitted plot, test quality, drainage condition, and stress convention.

Repository layout

MohrCoulomb-SimplexFit/
├── src/mohrcoulomb_simplexfit/  # numerical core, CLI, I/O, and plotting
├── tests/                       # analytical, numerical, and CLI verification
├── data/input/                  # local CSV entry point for interactive launchers
├── examples/                    # tracked reproducible datasets
├── docs/                        # method notes and README figures
├── scripts/                     # reproducible documentation-figure generation
├── run_fit.ps1                  # Windows interactive launcher
└── run_fit.sh                   # Bash interactive launcher

See the numerical method note, input-data guide, and contribution guide for further details.

Method references

1. Nelder, J. A., and Mead, R. (1965). “A simplex method for function minimization.” The Computer Journal, 7(4), 308–313. doi:10.1093/comjnl/7.4.308.
2. Labuz, J. F., and Zang, A. (2012). “Mohr–Coulomb Failure Criterion.” Rock Mechanics and Rock Engineering, 45, 975–979. doi:10.1007/s00603-012-0281-7.

If this software contributes to published work, report the version, input stress convention and unit, exact CSV data, fitted parameters, objective/RMSE, and convergence diagnostics. A formal software citation file can be added once project authorship and repository URL are finalized.

License

MohrCoulomb-SimplexFit is released under the MIT License.