Rbto tutorial

Project.toml

@@ -31,6 +31,7 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Preconditioners = "af69fa37-3177-5a40-98ee-561f696e4fcd"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+ReliabilityOptimization = "55eddd50-7f45-4398-96cf-6a37e2b16f80"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"

docs/Project.toml

@@ -1,11 +1,17 @@
 [deps]
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Nonconvex = "01bcebdf-4d21-426d-b5c4-6132c1619978"
+NonconvexMMA = "d3d89cbb-4ecd-4604-818d-8d1ff343e4da"
 NonconvexTOBS = "6c0b5230-d4c9-466e-bfd4-b31e6272ab65"
+ReliabilityOptimization = "55eddd50-7f45-4398-96cf-6a37e2b16f80"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 TopOpt = "53a1e1a5-51bb-58a9-8a02-02056cc81109"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

docs/make.jl

@@ -9,7 +9,15 @@ include("generate.jl")
 
 GENERATED_EXAMPLES = [
     joinpath("examples", f) for
-    f in ("simp.md", "beso.md", "geso.md", "csimp.md", "global_stress.md", "TOBS.md")
+    f in (
+        "simp.md",
+        "beso.md",
+        "geso.md", 
+        "csimp.md",
+        "global_stress.md",
+        "TOBS.md",
+        "RBTO.md",
+    )
 ]
 
 bib = CitationBibliography(joinpath(@__DIR__, "biblio", "ref.bib"))

docs/src/literate/RBTO.jl

@@ -0,0 +1,56 @@
+# # Multimaterial RBTO example 
+
+# [Reliability-based topology optimization](https://link.springer.com/article/10.1007/s00158-010-0518-6) (RBTO) allows the search for an optimal design while considering uncertainty in the problem, such as material properties. In the following, an example of RBTO with two materials is developed. First, packages and basic problem definitions are needed.
+
+# Packages and definitions.
+
+using ReliabilityOptimization, Test, ChainRulesCore
+using TopOpt, Zygote, FiniteDifferences, Nonconvex
+
+Nonconvex.@load MMA
+
+const v = 0.3 # Poisson’s ratio
+const densities = [0.0, 1.0] # for mass calc
+const f = 1.0 # downward force
+const nmats = 2
+const V = 0.4 # volume fraction
+const problemSize = (4, 4) # size of rectangular mesh
+const elSize = (1.0, 1.0); # size of QUAD4 elements
+# Point load cantilever problem to be solved
+problem = PointLoadCantilever(Val{:Linear}, problemSize, elSize, 1.0, v, f)
+ncells = TopOpt.getncells(problem) # Number of elements in mesh
+solver = FEASolver(Direct, problem; xmin = 0.0)
+filter = DensityFilter(solver; rmin = 3.0) # filter to avoid checkerboarding
+comp = Compliance(solver) # function that returns compliance
+penalty = TopOpt.PowerPenalty(3.0) # SIMP penalty
+const avgEs = [1e-6, 0.5]
+logEs = MvNormal(log.(avgEs), Matrix(Diagonal(abs.(log.(avgEs) * 0.1))))
+# 'Original' function. At least one input is random. In this example, Es is the random input.
+function compObj(y, logEs)
+    penalty1 = TopOpt.PowerPenalty(3.0)
+    interp1 = MaterialInterpolation(exp.(logEs), penalty1)
+    x = tounit(MultiMaterialVariables(y, nmats))
+    _E = interp1(filter(x))
+    dsd = comp(_E)
+    return dsd
+end
+penalty2 = TopOpt.PowerPenalty(1.0)
+interp2 = MaterialInterpolation(densities, penalty2)
+rf = RandomFunction(compObj, logEs, FORM(RIA())) # wrap original function in RandomFunction struct
+x0 = fill(V, ncells * (length(logEs) - 1)) # initial homogeneous distribution of pseudo-densities
+# call wrapper with example input. This returns a propability distribution of the objective for current point
+d = rf(x0)
+
+constr = y -> begin
+    _rhos = interp2(MultiMaterialVariables(y, nmats))
+    return sum(_rhos.x) / ncells - 0.4 # elements have unit volumes
+end
+function obj(x) # objective for TO problem
+    dist = rf(x)
+    mean(dist)[1] + 2 * sqrt(cov(dist)[1, 1])
+end
+
+m = Model(obj) # create optimization model
+addvar!(m, zeros(length(x0)), ones(length(x0))) # setup optimization variables
+Nonconvex.add_ineq_constraint!(m, constr) # setup volume inequality constraint
+@time r = Nonconvex.optimize(m, MMA87(), x0; options = MMAOptions())