Convex optimization-based constitutive update#
Elastoplastic return mapping via convex optimization#
For a specific class of materials, constitutive update can be also formulated as solving a convex optimization problem. In the case of associated plasticity, this amounts to projecting the elastic trial state onto the yield surface, [Bruno et al., 2020, Krabbenhøft et al., 2007]. Ignoring hardening for the sake of simplicity, such a projection can be formulated as follows:
where \(\bsig_\text{el} = \bsig_\text{old} + \mathbb{C}:\Delta\beps\) is the elastic predictor computed from the previous stress state \(\bsig_\text{old}\) and the current strain increment \(\Delta\beps\).
Writing the optimality conditions of (10) results in the system of nonlinear equations arising in the elastoplastic return mapping procedure. However, when considering non-smooth yield surfaces or yield conditions formed by the intersection of many yield surfaces, such local Newton procedures may encounter difficulties of convergence, due to the lack of differentiability of the yield surface.
In such cases, it may become interesting to consider using dedicated tools to solve the optimization problem (10) directly. At first sight, this problem is a quadratic optimization problem with nonlinear constraints. However, the yield function \(f(\bsig)\) can very often be reformulated using standard cones such as the positive orthant, the second-order Lorentz cone, the cone of positive-definite matrices, etc [Bisbos and Pardalos, 2007]. In this case, we can leverage more efficient convex optimization solvers such as primal-dual interior point solvers [Andersen et al., 2003].
CVXPY#
cvxpy is an open source Python-embedded modeling language for convex optimization problems, enabling to express such problems in a natural way using predefined atomic functions [Agrawal et al., 2018, Diamond and Boyd, 2016]. It offers a wide range of open-source and commercial solvers. Some of them are also capable of computing sensitivities of the optimization problem to derive the underlying tangent operator [Agrawal et al., 2019]. This is however not explored here at the moment.
See also
For more details, see also Andrey’s Latyshev work at a-latyshev/convex-plasticity.
Implementation#
We describe below the implementation of a base class called CvxPyMaterial. For simplicity, such materials are limited to plane stress elastoplasticity without hardening. We first initialize the class which inherits from the generic Material and enforce the plane dimension of stresses and strains.
from dolfinx_materials.material import Material
import cvxpy as cp
import numpy as np
class CvxPyMaterial(Material):
def __init__(self, elastic_model, **kwargs):
super().__init__()
self.elastic_model = elastic_model
# Handle any additional keyword arguments passed
for key, value in kwargs.items():
setattr(self, key, value)
self.C = np.array(self.elastic_model.compute_C_plane_stress())
self.set_cvxpy_model()
@property
def gradients(self):
return {"Strain": 3}
@property
def fluxes(self):
return {"Stress": 3}
The implementation of the optimization problem with cvxpy is done in the set_cvxpy_model method. The optimization variable is a stress state of dimension 3. The elastic stress is defined as a cvxpy.Parameter which are symbolic representations of constants. Using parameters allows to modify the values of constants without reconstructing the entire problem. However, this places some restriction on the way the problem is defined. These aspects are described in the Disciplined Parametrized Programming section of cvxpy documentation. Here, we define the objective function obj corresponding to \(\dfrac{1}{2}(\bsig-\bsig_\text{el}):\mathbb{C}^{-1}:(\bsig-\bsig_\text{el})\). The constraints are provided by the yield_constraints method which returns an empty list at the moment. Concrete implementation of subclasses of CvxPyMaterial should implementation the constraints corresponding to \(f(\bsig) \leq 0\).
def set_cvxpy_model(self):
self.sig = cp.Variable((3,))
self.sig_el = cp.Parameter((3,))
obj = 0.5 * cp.quad_form(self.sig - self.sig_el, np.linalg.inv(self.C))
self.prob = cp.Problem(cp.Minimize(obj), self.yield_constraints(self.sig))
def yield_constraints(self, Sig):
return []
Finally, the constitutive update rule is defined by computing the elastic predictor, affecting its value to the cvxpy.Parameter object and solving the problem. At the moment, this function returns the elastic operator only.
def constitutive_update(self, eps, state, dt):
eps_old = state["Strain"]
deps = eps - eps_old
sig_old = state["Stress"]
sig_pred = sig_old + self.C @ deps
self.sig_el.value = sig_pred
self.prob.solve()
state["Strain"] = eps
state["Stress"] = self.sig.value
return self.C, state
References#
Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. Advances in neural information processing systems, 2019.
Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd. A rewriting system for convex optimization problems. Journal of Control and Decision, 5(1):42–60, 2018.
Erling D Andersen, Cornelis Roos, and Tamas Terlaky. On implementing a primal-dual interior-point method for conic quadratic optimization. Mathematical Programming, 95(2):249–277, 2003.
Christos D Bisbos and Panos M Pardalos. Second-order cone and semidefinite representations of material failure criteria. Journal of Optimization Theory and Applications, 134(2):275–301, 2007.
Hugo Bruno, Guilherme Barros, Ivan FM Menezes, and Luiz Fernando Martha. Return-mapping algorithms for associative isotropic hardening plasticity using conic optimization. Applied Mathematical Modelling, 78:724–748, 2020.
Steven Diamond and Stephen Boyd. CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1–5, 2016.
K Krabbenhøft, AV Lyamin, and SW Sloan. Formulation and solution of some plasticity problems as conic programs. International Journal of Solids and Structures, 44(5):1533–1549, 2007.