１０本トラスの最適化問題を差分法を用いて求めたいです。アルゴリズム通りに計算をしてはずなのですが、答えが合いません。コードを見て頂けませんでしょうか。

前提・実現したいこと

１０本のトラスの最適化問題を解いております。
有限差分法を用いて、計算をし、scipy.minimize にて、最適化をしています。

発生している問題

アルゴリズム通りにコードを書いたのですが、理想の答え（８３０）が出ませんでした。
現在は、232が結果として表示されます。
考えられる間違いとしては、
・スケーリングの値
・有限差分法の計算の間違い
・勾配の引数の値の間違い
があるのではないか考えたのですが、何度確認してみてもコードに間違いはない、と思ってしまい、
どこで間違えたのか、理解ができません。
教えて頂けないでしょうか。

該当のソースコード

1つ目は、目的関数の設定と、有限差分法を用いて、勾配を求めるコードです。

python
1import numpy as np
2#from truss_visualizer import cool_plot
3from copy import deepcopy
4from math import sin, cos, pi
5
6#ーーーーーーーー有限差分法ーーーーーーーーーーーー
7def fd(A):
8    dmass_dA = np.empty_like(A)
9    dstress_dA = np.empty((A.shape[0],A.shape[0]))
10    h=1e-4
11    for i in range(A.shape[0]):
12        Xp = deepcopy(A)
13        Xp[i] += h
14        mass_p, stress_p = trussf(Xp)
15
16        Xm = deepcopy(A)
17        Xm[i] -= h
18        mass_m, stress_m = trussf(Xm)
19
20        dmass_dA[i] = (mass_p - mass_m) / (2.0 * h)
21        dstress_dA[i] = (stress_p - stress_m) / (2.0 * h)
22
23    return dmass_dA, dstress_dA
24
25#ーーーーーーーー目的関数の計算（初期値がコードに含まれています）ーーーーーーーーーーーー
26def trussf(A):
27
28    # --- setup 10 bar truss ----
29    n=10
30    nbar = np.size(A)
31    E=70e9* np.ones(n)
32    rho=2720 * np.ones(n)
33    P=5e5
34    Len_side=10
35    Len_diagonal=10*np.sqrt(2)
36    L=np.concatenate(([Len_side] * 6, [Len_diagonal] * 4), axis=0)
37    start=np.array([5, 3, 6, 4, 4, 2, 4, 6, 2, 4])
38    finish=np.array([3, 1, 4, 2, 3, 1, 5, 3, 3, 1])
39    phi = np.array([0, 0, 0, 0, pi/2, pi/2, 3*pi/4, pi/4, 3*pi/4, pi/4])
40    rigid = [False, False, False, False, True, True]
41    Fx= np.zeros(6, dtype=float)
42    Fy= np.array([0.0, -P, 0.0, -P, 0.0, 0.0])
43
44    # --- call truss function ----
45    mass, stress = truss(start, finish, phi, A, L, E, rho, Fx, Fy, rigid)
46
47
48    return mass, stress
49
50def bar(E, A, L, phi):
51    """Computes the stiffness and stress matrix for one element
52
53    Parameters
54    ----------
55    E : float
56        modulus of elasticity
57    A : float
58        cross-sectional area
59    L : float
60        length of element
61    phi : float
62        orientation of element
63
64    Outputs
65    -------
66    K : 4 x 4 ndarray
67        stiffness matrix
68    S : 1 x 4 ndarray
69        stress matrix
70
71    """
72
73    # rename
74    c = cos(phi)
75    s = sin(phi)
76
77    # stiffness matrix
78    k0 = np.array([[c**2, c*s], [c*s, s**2]])
79    k1 = np.hstack([k0, -k0])
80    K = E*A/L*np.vstack([k1, -k1])
81
82    # stress matrix
83    S = E/L*np.array([[-c, -s, c, s]])
84
85    return K, S
86
87def node2idx(node, DOF):
88    """Computes the appropriate indices in the global matrix for
89    the corresponding node numbers.  You pass in the number of the node
90    (either as a scalar or an array of locations), and the degrees of
91    freedom per node and it returns the corresponding indices in
92    the global matrices
93
94    """
95
96    idx = np.array([], dtype=np.int)
97
98    for i in range(len(node)):
99
100        n = node[i]
101        start = DOF*(n-1)
102        finish = DOF*n
103
104        idx = np.concatenate((idx, np.arange(start, finish, dtype=np.int)))
105
106    return idx
107
108#ーーーーーーーー目的関数の計算（初期値がコードに含まれていません）ーーーーーーーーーーーー
109def truss(start, finish, phi, A, L, E, rho, Fx, Fy, rigid):
110    """Computes mass and stress for an arbitrary truss structure
111
112    Parameters
113    ----------
114    start : ndarray of length nbar
115        index of start of bar (1-based indexing) start and finish can be in any order as long as consistent with phi
116    finish : ndarray of length nbar
117        index of other end of bar (1-based indexing)
118    phi : ndarray o f length nbar (radians)
119        defines orientation or bar
120    A : ndarray of length nbar
121        cross-sectional areas of each bar
122    L : ndarray of length nbar
123        length of each bar
124    E : ndarray of length nbar
125        modulus of elasticity of each bar
126    rho : ndarray of length nbar
127        material density of each bar
128    Fx : ndarray of length nnode
129        force in the x-direction at each node
130    Fy : ndarray of length nnode
131        force in the y-direction at each node
132    rigid : list(boolean) of length nnode
133        True if node_i is rigidly constrained
134
135    Outputs
136    -------
137    mass : float
138        mass of the entire structure
139    stress : ndarray of length nbar
140        stress of each bar
141
142    """
143
144    n = len(Fx)  # number of nodes
145    DOF = 2  # number of degrees of freedom
146    nbar = len(A)  # number of bars
147
148    # mass
149    mass = np.sum(rho*A*L)
150
151    # stiffness and stress matrices
152    K = np.zeros((DOF*n, DOF*n), dtype=A.dtype)
153    S = np.zeros((nbar, DOF*n), dtype=A.dtype)
154
155    for i in range(nbar):  # loop through each bar
156
157        # compute submatrix for each element
158        Ksub, Ssub = bar(E[i], A[i], L[i], phi[i])
159
160        # insert submatrix into global matrix
161        idx = node2idx([start[i], finish[i]], DOF)  # pass in the starting and ending node number for this element
162        K[np.ix_(idx, idx)] += Ksub
163        S[i, idx] = Ssub
164
165    # applied loads
166    F = np.zeros((n*DOF, 1))
167
168    for i in range(n):
169        idx = node2idx([i+1], DOF)  # add 1 b.c. made indexing 1-based for convenience
170        F[idx[0]] = Fx[i]
171        F[idx[1]] = Fy[i]
172
173
174    # boundary condition
175    idx = np.squeeze(np.where(rigid))
176    remove = node2idx(idx+1, DOF)  # add 1 b.c. made indexing 1-based for convenience
177
178    K = np.delete(K, remove, axis=0)
179    K = np.delete(K, remove, axis=1)
180    F = np.delete(F, remove, axis=0)
181    S = np.delete(S, remove, axis=1)
182
183    # solve for deflections
184    d = np.linalg.solve(K, F)
185
186    # compute stress
187    stress = np.dot(S, d).reshape(nbar)
188
189    return mass, stress
190
191#ーーーーーーーー最初に呼び出される目的関数ーーーーーーーーーーー
192def tenbartruss(A, grad_type):
193    """This is the subroutine for the 10-bar truss.  You will need to complete it.
194
195    Parameters
196    ----------
197    A : ndarray of length 10
198        cross-sectional areas of all the bars
199    grad_type : string (optional)
200        gradient type.  'FD' for finite difference, 'CS' for complex step,
201        'AJ' for adjoint
202
203    Outputs
204    -------
205    mass : float
206        mass of the entire structure
207    stress : ndarray of length 10
208        stress of each bar
209    dmass_dA : ndarray of length 10
210        derivative of mass w.r.t. each A
211    dstress_dA : 10 x 10 ndarray
212        dstress_dA[i, j] is derivative of stress[i] w.r.t. A[j]
213
214    """
215
216    # --- setup 10 bar truss ----
217    n=10
218    nbar = np.size(A)
219    E=70e9* np.ones(n)
220    rho=2720 * np.ones(n)
221    P=5e5
222    Len_side=10
223    Len_diagonal=10*np.sqrt(2)
224    L=np.concatenate(([Len_side] * 6, [Len_diagonal] * 4), axis=0)
225    start=np.array([5, 3, 6, 4, 4, 2, 4, 6, 2, 4])
226    finish=np.array([3, 1, 4, 2, 3, 1, 5, 3, 3, 1])
227    phi = np.array([0, 0, 0, 0, pi/2, pi/2, 3*pi/4, pi/4, 3*pi/4, pi/4])
228    rigid = [False, False, False, False, True, True]
229    Fx= np.zeros(6, dtype=float)
230    Fy= np.array([0.0, -P, 0.0, -P, 0.0, 0.0])
231
232    # --- call truss function ----
233    mass, stress = truss(start, finish, phi, A, L, E, rho, Fx, Fy, rigid)
234
235    # --- compute derivatives for provided grad_type ----
236
237    if grad_type == 'FD':
238        dmass_dA, dstress_dA = fd(A)
239
240    return mass*1e-6, stress*1e-6, dmass_dA*1e-6, dstress_dA*1e-6
241

（呼び出すコードです）2つ目は、1つ目で求めた勾配や目的関数の値を、scipyを用いて最適化するコードです。

python
1
2import numpy as np
3from scipy.optimize import minimize
4from scipy.optimize import NonlinearConstraint
5from truss import tenbartruss
6import sys
7from IPython.core.debugger import set_trace
8from truss_FD import truss
9
10grad_method='FD'
11num_evals = 0
12mass = 0
13stress = 0
14all_functions = []
15max_const_vio = []
16
17def solve(x):
18    global num_evals
19    num_evals += 1
20    mass, stress, dmdA, dsdA = tenbartruss(x,grad_method)
21    #mass, stress, dmdA, dsdA = truss(x)
22    return mass
23
24
25def constraintAll(x):
26    global num_evals
27    num_evals += 1
28    mass, stress, dmdA, dsdA= tenbartruss(x,grad_method)
29    #mass, stress, dmdA, dsdA = truss(x)
30
31    stress_ub = 170e6*1e-6 - stress
32    stress_ub[8] = 520e6*1e-6 - stress[8]
33
34    stress_lb = stress + 170e6*1e-6
35    stress_lb[8] = stress[8] + 520e6*1e-6
36
37
38    return np.hstack((stress_ub,stress_lb))
39
40def get_dsdA(x):
41    global num_evals
42    num_evals += 1
43    mass, stress, dmdA, dsdA = tenbartruss(x,grad_method)
44    #mass, stress, dmdA, dsdA = truss(x)
45    return np.vstack((-dsdA, dsdA))
46
47
48def get_dmdA(x):
49    global num_evals
50    num_evals += 1
51    mass, stress, dmdA, dsdA = tenbartruss(x,grad_method)
52    #mass, stress, dmdA, dsdA = truss(x)
53    return dmdA
54
55def callbackF2(Xi):
56    temp = 1
57    global all_functions, max_const_vio
58    print('ok',solve(Xi))
59    all_functions.append(solve(Xi))
60    max = np.max(np.abs(constraintAll(Xi)))
61
62    max_const_vio.append(max)
63
64
65all_constraints = {'type':'ineq', 'fun':constraintAll, 'jac':get_dsdA}
66callback = callbackF2
67
68#### Solve
69start = np.ones(10) * 0.1
70boundary=((5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None),(5e-5,None))
71fit = minimize(solve,
72    start,
73    method="SLSQP",
74    jac=get_dmdA,
75    bounds=boundary,
76    constraints=all_constraints,
77    tol=1e-12,
78    options={'disp':True,'maxiter':500},
79    callback=callback)
80
81#### Output values for debugging
82print("Cross Sectional Area")
83print(fit.x)
84print("-----------------------")
85print("minima and stress")
86print(tenbartruss(fit.x,grad_method)
87print("num of iteration")
88print(num_evals)
89print("-----------------------")
90
91

試したこと

1つ目のコードの最後の行のスケーリングの値を少しずつ調整しました。
コードを１ステップずつ実行し、どこで計算が違って居るか、確認しようとしましたが、見つかりませんでした。
アルゴリズムの計算上では、８３０と結果が表示されるはずなのですが、スケーリングを変えても、差分法のコードを変えても、上の値は出せませんでした。
私のプログラムのおかしい点を、教えて頂けないでしょうか。