diff --git a/Project.toml b/Project.toml
index 97010245..a2f06b37 100644
--- a/Project.toml
+++ b/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"
diff --git a/docs/Project.toml b/docs/Project.toml
index 8a65511c..ad1faee1 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/docs/make.jl b/docs/make.jl
index e0be28b8..e46d1603 100644
--- a/docs/make.jl
+++ b/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"))
diff --git a/docs/src/literate/RBTO.jl b/docs/src/literate/RBTO.jl
new file mode 100644
index 00000000..55262f61
--- /dev/null
+++ b/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())
\ No newline at end of file