In [1]:
%matplotlib inline
execfile('../../shim_preamble.py')
In [2]:
cd ~/Workspace/gm2/studies/field-shimming/field-modeling
Optimizing the Field with Shim Adjustments v0.6¶
This is an improvement to the last iteration.
Changes from v0.5¶
- lets edge shims change assymetrically
- ability to load/upload shim adjustments to google spreadsheet
- print/store plans in mils for top hats and turns for wedges
- refactored 'tophat' in favor of 'top_hat' everywhere
- redundant implementation code is now loaded from a single file for all notebooks
Shim Models¶
In [3]:
# Shim models for The optimization routine.
def current_model(x, pos, mp_id):
"""Current model creates a constant offset for the average field."""
if mp_id == 0:
amp = 1.0
else:
amp = 0.0
return amp * pos * np.ones(x.shape)
def top_hat_model(x, pos, yoke_id, top_hat_id, mp_id):
"""Define model for top hat shims."""
phi_0 = yoke_id * 30.0 - 7.5 + 15.0 * top_hat_id
phi = (x - phi_0 + 180.0) % 360 - 180.0
if mp_id == 0:
amp = -372.0 # -186.0 for only uppers
else:
amp = 0.0
return amp * pos * np.exp(-(phi / 20.74)**2.0)
def wedge_model(x, pos, pole_id, wedge_id, mp_id):
"""Define model for wedge shims."""
phi_0 = 10.0 * pole_id - 5.0 + 10.0 / 12.0 * (wedge_id + 0.5)
phi = (x - phi_0 + 180.0) % 360 - 180.0
if mp_id == 0:
amp = 4.0
elif mp_id == 1:
amp = -0.11
elif mp_id == 2:
amp = 0.0 # -0.28 only if not symmetric
else:
amp = 0.0
return amp * pos * np.exp(-(phi / 5.139)**2.0)
def inner_upper_edge_model(x, pos, pole_id, edge_id, mp_id):
"""Define model for inner upper edge shims."""
phi_0 = 10.0 * pole_id
y = np.zeros(x.shape)
for x0 in np.linspace(phi_0 - 5.0, phi_0 + 5.0, 20):
phi = (x - x0 + 180.0) % 360 - 180.0
y += np.exp(-(phi / 0.5)**2.0)
y /= y.max()
if mp_id == 1:
amp = -38.0
elif mp_id == 2:
amp = -216.0
elif mp_id == 3:
amp = -31.2
elif mp_id == 4:
amp = -81.2
elif mp_id == 5:
amp = -16.0
elif mp_id == 6:
amp = -19.6
else:
amp = 0.0
return amp * pos * y
def inner_lower_edge_model(x, pos, pole_id, edge_id, mp_id):
"""Define model for inner lower edge shims."""
phi_0 = 10.0 * pole_id
y = np.zeros(x.shape)
for x0 in np.linspace(phi_0 - 5.0, phi_0 + 5.0, 20):
phi = (x - x0 + 180.0) % 360 - 180.0
y += np.exp(-(phi / 0.5)**2.0)
y /= y.max()
if mp_id == 1:
amp = -38.0
elif mp_id == 2:
amp = +216.0
elif mp_id == 3:
amp = -31.2
elif mp_id == 4:
amp = +81.2
elif mp_id == 5:
amp = -16.0
elif mp_id == 6:
amp = +19.6
else:
amp = 0.0
return amp * pos * y
def outer_edge_model(x, pos, pole_id, edge_id, mp_id):
"""Define model for outer edge shims."""
phi_0 = 10.0 * pole_id
y = np.zeros(x.shape)
for x0 in np.linspace(phi_0 - 5.0, phi_0 + 5.0, 20):
phi = (x - x0 + 180.0) % 360 - 180.0
y += np.exp(-(phi / 0.5)**2.0)
y /= y.max()
if mp_id == 1:
amp = 33.2
elif mp_id == 2:
amp = -195.2
elif mp_id == 3:
amp = -26.4 * 2
elif mp_id == 4:
amp = 74.0 * 0
elif mp_id == 5:
amp = -13.2 * 2
elif mp_id == 6:
amp = -17.2 * 0
else:
amp = 0.0
return amp * pos * y
Load In¶
We want to load the field data to be optimized, and set up the shim models with initial conditions.
In [4]:
datafile = 'data/extracted/full_scan_034.root'
if datafile == 'data/extracted/full_scan_028.root':
phi_nmr_offset = 0.0
elif datafile == 'data/extracted/full_scan_032.root':
phi_nmr_offset = 0.1
elif datafile == 'data/extracted/full_scan_033.root':
phi_nmr_offset = 0.36
elif datafile == 'data/extracted/full_scan_034.root':
phi_nmr_offset = 0.36
else:
phi_nmr_offset = 1.42
npoints = 2000
max_iter = 1000
wedge_turns_per_mm = 1.6
upload_shim_adjustments = False
enable_shim = {}
enable_shim['current'] = True
enable_shim['top_hat'] = True
enable_shim['wedge'] = True
enable_shim['edge'] = False
# Define the current top_hat positions
# top_hat_pos = [0.0, 0.0, # A1/A2
# 0.0, 0.0, # B1/B2
# 1.0, 1.0, # C1/C2
# 0.0, 0.5, # D1/D2
# 0.0, 0.0, # E1/E2
# 1.0, 1.0, # F1/F2
# 0.0, 0.0, # G1/G2
# 0.0, 1.5, # H1/H2
# 1.0, 1.0, # I1/I2
# 1.0, 0.5, # J1/J2
# 0.0, 0.0, # K1/K2
# 0.0, 1.0] # L1/L2
top_hat_pos = get_current_top_hat_settings()
wedge_pos = get_current_wedge_settings()
current_lb = 49
current_ub = 50
top_hat_lb = 0.0
top_hat_ub = 2.0
wedge_lb = -10.0
wedge_ub = 10.0
edge_lb = 0.0
edge_ub = 0.2
# Constants
num_yokes = 12
num_poles = 72
num_top_hats = 24
num_edges = num_poles * 2
num_wedges = 864 / 2 # Move top/bottom symmetrically
# Set the relative weights for removing multipole
multipole_wt = np.array([1.0, 2.5, 2.5, 10.0, 10.0, 10.0, 10.0])
num_mp = len(multipole_wt)
mp_name = {}
mp_name[0] = 'Dipole'
mp_name[1] = 'Normal Quadrupole'
mp_name[2] = 'Skew Quadrupole'
mp_name[3] = 'Normal Sextupole'
mp_name[4] = 'Skew Sextupole'
mp_name[5] = 'Normal Octupole'
mp_name[6] = 'Skew Octupole'
In [5]:
# Load the data.
f = rt.TFile(datafile)
t = f.Get('t')
fs_phi = np.empty(t.GetEntries())
fs_bfield = np.empty([num_mp, t.GetEntries()])
for i in xrange(t.GetEntries()):
t.GetEntry(i)
fs_phi[i] = (t.phi_2 - phi_nmr_offset) % 360.0
for j in xrange(num_mp):
fs_bfield[j, i] = t.multipole[j]
phi = np.linspace(0.0, 360.0, npoints)
bfield = np.empty([num_mp, npoints])
for i in xrange(num_mp):
bfield[i] = griddata(np.hstack((fs_phi[-100:] - 360.0, fs_phi, fs_phi[:100] + 360.0)),
np.hstack((fs_bfield[i, -100:], fs_bfield[i], fs_bfield[i, :100])), phi)
# Allocated the the shim matrix.
num_knobs = 1 + num_top_hats + num_wedges + 2 * num_edges
A = np.empty([num_knobs, num_mp, npoints])
lb = np.empty(A.shape[0])
ub = np.empty(A.shape[0])
# Build the shim matrix, constant term first (magnet current)
if enable_shim['current']:
for i in xrange(num_mp):
A[0, i] = current_model(phi, 1.0, i)
lb[0] = current_lb
ub[0] = current_ub
else:
print "Loading zeros for current."
for i in xrange(num_mp):
A[0, i] = np.zeros(A[0, i].shape)
lb[0] = -0.00001
ub[0] = +0.00001
idx = 0
# Next input templates for the top_hat.
if enable_shim['top_hat']:
for i in xrange(num_yokes):
for j in xrange(num_top_hats / num_yokes):
idx += 1
lb[idx] = top_hat_lb - top_hat_pos[2 * i + j] # allows for [0mm, 2mm]
ub[idx] = top_hat_ub - top_hat_pos[2 * i + j]
for k in xrange(num_mp):
A[idx, k] = multipole_wt[k] * top_hat_model(phi, 1.0, i, j, k)
else:
print "Loading zeros for top_hat."
for i in xrange(num_yokes):
for j in xrange(num_top_hats / num_yokes):
idx += 1
lb[idx] = -0.00001
ub[idx] = +0.00001
for k in xrange(num_mp):
A[idx, k] = np.zeros(A[idx, k].shape)
# Now add the wedge shims.
if enable_shim['wedge']:
for i in xrange(num_poles / 2):
for j in xrange(num_wedges / (num_poles / 2)):
idx += 1
lb[idx] = wedge_lb - wedge_pos[i * 12 + j]
ub[idx] = wedge_ub - wedge_pos[i * 12 + j]
for k in xrange(num_mp):
A[idx, k] = multipole_wt[k] * wedge_model(phi, 1.0, i, j, k)
else:
print "Loading zeros for wedges."
for i in xrange(num_poles / 2):
for j in xrange(num_wedges / (num_poles / 2)):
idx += 1
lb[idx] = -0.00001
ub[idx] = +0.00001
for k in xrange(num_mp):
A[idx, k] = np.zeros(A[idx, k].shape)
if enable_shim['edge']:
for i in xrange(num_poles / 2):
for j in xrange(num_edges / (num_poles / 2)):
idx += 1
lb[idx] = edge_lb
ub[idx] = edge_ub
for k in xrange(num_mp):
A[idx, k] = multipole_wt[k] * inner_edge_model(phi, 1.0, i, j, k)
for i in xrange(num_poles / 2):
for j in xrange(num_edges / (num_poles / 2)):
idx += 1
lb[idx] = edge_lb
ub[idx] = edge_ub
for k in xrange(num_mp):
A[idx, k] = multipole_wt[k] * outer_edge_model(phi, 1.0, i, j, k)
else:
print "Loading zeros for edges."
for i in xrange(num_poles / 2):
for j in xrange(num_edges / (num_poles / 2)):
idx += 1
lb[idx] = -0.00001
ub[idx] = +0.00001
for k in xrange(num_mp):
A[idx, k] = np.zeros(A[idx, k].shape)
for i in xrange(num_poles / 2):
for j in xrange(num_edges / (num_poles / 2)):
idx += 1
lb[idx] = -0.00001
ub[idx] = +0.00001
for k in xrange(num_mp):
A[idx, k] = np.zeros(A[idx, k].shape)
A = A.reshape([num_knobs, npoints * num_mp])
# Set up the field.
B0 = np.empty([num_mp, npoints])
field = np.empty([num_mp, npoints])
for k in xrange(num_mp):
field[k] = multipole_wt[k] * bfield[k]
if k == 0:
B0[k, :] = bfield[k].mean() * np.ones(bfield[k].shape)
else:
B0[k, :] = np.zeros(bfield[k].shape)
B0 = B0.flatten()
field = field.flatten()
Optimization¶
Actually perform the least square minimization.
In [6]:
res = lsq_linear(A.T, B0 - field, (lb, ub), method='trf', max_iter=max_iter, verbose=1)
x = res.x
In [7]:
# ff1 = Ridge(alpha=0.1, max_iter=5000, normalize=True)
# ff1.fit(A.T, B0 - field)
# x = ff1.coef_
In [8]:
# ff2 = Lasso(alpha=0.1, max_iter=5000)
# ff2.fit(A.T, B0 - field)
# x = ff1.coef_
In [9]:
print "dB = %.3f" % x[0]
y = np.dot(A.T, x) + field
y_bfield = y.reshape([num_mp, npoints])
y_shims = np.empty([3, num_mp, npoints])
indices = np.array(range(1, 1 + num_top_hats))
y_top_hats = np.dot(A[indices, :].T, x[indices])
y_top_hats = y_top_hats.reshape([num_mp, npoints])
indices = np.array(range(1 + num_top_hats, 1 + num_top_hats + num_wedges))
y_wedges = np.dot(A[indices, :].T, x[indices])
y_wedges = y_wedges.reshape([num_mp, npoints])
indices = np.array(range(1 + num_top_hats + num_wedges, 1 + num_top_hats + num_wedges + num_edges * 2))
y_edges = np.dot(A[indices, :].T, x[indices])
y_edges = y_edges.reshape([num_mp, npoints])
for i in xrange(num_mp):
if i == 0:
y_bfield[i] -= x[0]
y_bfield[i] /= multipole_wt[i]
y_top_hats[i] /= multipole_wt[i]
y_wedges[i] /= multipole_wt[i]
y_edges[i] /= multipole_wt[i]
Resulting Field¶
In [10]:
for i in xrange(num_mp):
print "Current %s avg: %.3f, std: %.3f" % (mp_name[i], bfield[i].mean(), bfield[i].std())
print "Optimal %s avg: %.3f, std: %.3f" % (mp_name[i], y_bfield[i].mean(), y_bfield[i].std())
# Draw the multipole with target lines and original data.
plt.scatter(phi, bfield[i], color=colors[0])
plt.scatter(phi, y_bfield[i], color=colors[2])
if i == 0:
draw_targets(y_bfield[i].mean(), 25.0)
else:
draw_targets(0.0, 10.0)
plt.ylim([-2 * y_bfield[i].std(), 2 * y_bfield[i].std()])
finish_plot(ylabel='%s [ppm]' % mp_name[i])
# Draw the multipole with target lines alone.
plt.scatter(phi, y_bfield[i], color=colors[2])
if i == 0:
draw_targets(y_bfield[i].mean(), 25.0)
else:
draw_targets(0.0, 10.0)
plt.ylim([-40, 40])
finish_plot(ylabel='%s [ppm]' % mp_name[i])
Field From Shims¶
I'm plotting the field change affected by the optimization model. The top hats are in red, the wedges are in green and the eges are in blue.
In [11]:
for i in xrange(num_mp):
plt.scatter(phi, y_top_hats[i], color=colors[0])
plt.scatter(phi, y_wedges[i], color=colors[2])
plt.scatter(phi, y_edges[i], color=colors[6])
finish_plot(ylabel='%s [ppm]' % mp_name[i])
Shim Positions¶
In [12]:
xtemp = x[1:num_top_hats + 1]
y1 = np.empty(xtemp.shape[0] * 100)
for i in xrange(y1.shape[0]):
y1[i] = xtemp[i / 100] + top_hat_pos[i / 100]
x1 = np.linspace(-15.0, 345.0, num_top_hats * 100) % 360.0
indices = x1.argsort()
plt.fill_between(x1[indices], y1[indices], color=colors[0])
finish_plot('top_hat shim position [mm]')
xtemp = x[1 + num_top_hats:1 + num_top_hats + num_wedges]
y2 = np.empty(xtemp.shape[0] * 100)
for i in xrange(y2.shape[0]):
y2[i] = xtemp[i / 100] + wedge_pos[i / 100]
x2 = np.linspace(-5, 355, num_wedges * 100) % 360.0
indices = x2.argsort()
plt.fill_between(x2[indices], y2[indices], color=colors[2])
finish_plot('wedge shim position [mm]')
plt.hist(x[1:num_top_hats + 1] + np.array(top_hat_pos), bins=40)
plt.xlabel('top_hat shim position [mm]')
plt.show()
plt.hist(x[1 + num_top_hats:1 + num_top_hats + num_wedges] + np.array(wedge_pos[:432]), bins=80)
plt.xlabel('wedge shim position [mm]')
plt.show()
mpl.rcParams['figure.figsize'] = (18, 3)
xtemp = x[1:num_top_hats + 1]
y1 = np.empty(xtemp.shape[0] * 100)
for i in xrange(y1.shape[0]):
y1[i] = xtemp[i / 100] + top_hat_pos[i / 100]
x1 = np.linspace(-15.0, 345.0, num_top_hats * 100) % 360.0
indices = x1.argsort()
plt.fill_between(x1[indices], y1[indices], color=colors[0])
finish_plot('top_hat shim position [mm]')
xtemp = x[1 + num_top_hats:1 + num_top_hats + num_wedges]
y2 = np.empty(xtemp.shape[0] * 100)
for i in xrange(y2.shape[0]):
y2[i] = xtemp[i / 100] + wedge_pos[i / 100]
x2 = np.linspace(-5.0, 355.0, num_wedges * 100) % 360.0
indices = x2.argsort()
plt.fill_between(x2[indices], y2[indices], color=colors[2])
finish_plot('wedge shim position [mm]')
mpl.rcParams['figure.figsize'] = (18, 6)
Shim Adjustment Plan¶
Based on the field optimization calculation the following changes to the shim positions are recommended.
In [13]:
display(Markdown('### Top Hat Changes'))
output = ['| Yoke %c | Top Hat 1 | Top Hat 2 |']
output.append('|-------------|------|------|')
output.append('| Change [mil] | %.2f | %.2f |')
top_hat_delta = (x[1:1 + num_top_hats]).reshape([num_top_hats / 2, 2])
if upload_shim_adjustments:
update_top_hat_settings(top_hat_delta)
for i, d in enumerate(top_hat_delta / 0.0254):
display(Markdown('\n'.join(output) % tuple([ord('A') + i, d[0], d[1]])))
display(Markdown('### Top Hat Position'))
output = ['| Yoke %c | Top Hat 1 | Top Hat 2 |']
output.append('|-------------|------|------|')
output.append('| Height [mil] | %.2f | %.2f |')
top_hat_abs = (x[1:1 + num_top_hats] + np.array(top_hat_pos)).reshape([num_top_hats / 2, 2])
for i, d in enumerate(top_hat_abs / 0.0254):
display(Markdown('\n'.join(output) % tuple([ord('A') + i, d[0], d[1]])))
display(Markdown('### Wedge Shim Adjustments'))
output = ['| Pole %i | W-01 | W-02 | W-03 | W-04 | W-05 | W-06 | W-07 | W-08 | W-09 | W-10 | W-11 | W-12 |']
output.append('|----------|------|------|------|------|------|------|------|------|------|------|------|------|')
output.append('| Pos [turns] | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f |')
wedge_delta = (x[1 + num_top_hats:1 + num_wedges + num_top_hats]).reshape([num_wedges / 12, 12])
if upload_shim_adjustments:
update_wedge_settings(np.hstack((wedge_delta, wedge_delta)))
for i, d in enumerate(wedge_delta / wedge_turns_per_mm):
display(Markdown('\n'.join(output) % tuple([i + 1, d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8], d[9], d[10], d[11]])))
display(Markdown('### Wedge Shim Positions'))
output = ['| Pole %i | W-01 | W-02 | W-03 | W-04 | W-05 | W-06 | W-07 | W-08 | W-09 | W-10 | W-11 | W-12 |']
output.append('|----------|------|------|------|------|------|------|------|------|------|------|------|------|')
output.append('| Pos [turns] | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f | %.2f |')
wedge_abs = (x[1 + num_top_hats:1 + num_wedges + num_top_hats] + np.array(wedge_pos[:432])).reshape([num_wedges / 12, 12])
for i, d in enumerate(wedge_abs / wedge_turns_per_mm):
display(Markdown('\n'.join(output) % tuple([i + 1, d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8], d[9], d[10], d[11]])))
In [ ]: