1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
#! /usr/bin/env python
import openturns as ot
import openturns.testing as ott
import sys
ot.TESTPREAMBLE()
ot.PlatformInfo.SetNumericalPrecision(3)
def progress(percent):
sys.stderr.write("-- progress=" + str(percent) + "%\n")
def stop():
sys.stderr.write("-- stop?\n")
return False
dim = 2
f = ot.SymbolicFunction(["x1", "x2"], ["10-5*(x1-3)^2-7*(x2-2)^2+0.1*(x1+x2)"])
startingPoint = [0.0] * dim
bounds = ot.Interval([-6.0, -6.0], [1.0, 7.0])
algoNames = ot.NLopt.GetAlgorithmNames()
for algoName in algoNames:
algo = ot.NLopt(algoName)
if algoName == "GN_ISRES":
algo.setMaximumConstraintError(1e-2)
for minimization in [False, True]:
for inequality in [False, True]:
for equality in [False, True]:
for bound in [False, True]:
if minimization and not bound:
continue
print(
"algo=",
algoName,
"minimization=",
minimization,
"bounds=",
bound,
"inequality=",
inequality,
"equality=",
equality,
)
problem = ot.OptimizationProblem(f)
problem.setMinimization(minimization)
if inequality:
# x1 <= 2
problem.setInequalityConstraint(
ot.SymbolicFunction(["x1", "x2"], ["2-x1"])
)
if equality:
# x2 = 4
problem.setEqualityConstraint(
ot.SymbolicFunction(["x1", "x2"], ["x2-4"])
)
if bound:
problem.setBounds(bounds)
try:
algo.setProblem(problem)
except Exception:
print("-- Not supported")
continue
algo.setMaximumCallsNumber(1000)
algo.setStartingPoint(startingPoint)
try:
algo.run()
except Exception as e:
print("-- ", e)
continue
result = algo.getResult()
x = result.getOptimalPoint()
print("x^=", x, "y^=", result.getOptimalValue())
if not inequality and not equality:
if not minimization:
if not bound:
x_ref = [3.0, 2.0]
else:
x_ref = [1.0, 2.0]
else:
x_ref = [-6.0, -6.0]
# these go to (-6,7)
if "NELDERMEAD" in algoName or "SBPLX" in algoName:
continue
ott.assert_almost_equal(x, x_ref, 4e-1, 1e-2)
elif equality:
ott.assert_almost_equal(x[1], 4.0, 4e-1, 1e-2)
elif inequality:
assert x[0] < 2.01, "!x1<=2.0"
# FORM
f = ot.SymbolicFunction(["E", "F", "L", "I"], ["-F*L^3/(3*E*I)"])
dim = f.getInputDimension()
mean = [50.0, 1.0, 10.0, 5.0]
sigma = ot.Point(dim, 1.0)
R = ot.IdentityMatrix(dim)
distribution = ot.Normal(mean, sigma, R)
vect = ot.RandomVector(distribution)
output = ot.CompositeRandomVector(f, vect)
event = ot.ThresholdEvent(output, ot.Less(), -3.0)
solver = ot.NLopt("LD_AUGLAG")
solver.setMaximumCallsNumber(400)
solver.setMaximumAbsoluteError(1.0e-10)
solver.setMaximumRelativeError(1.0e-10)
solver.setMaximumResidualError(1.0e-10)
solver.setMaximumConstraintError(1.0e-10)
algo = ot.FORM(solver, event, mean)
algo.run()
result = algo.getResult()
beta = result.getGeneralisedReliabilityIndex()
print("beta=%.6f" % beta)
ott.assert_almost_equal(beta, 1.009255)
|
