VIEWS: 6 PAGES: 71 POSTED ON: 9/29/2011 Public Domain
Patrick Gloster Summer Project 2008: Compact Straightness Monitor Simulation/Multilateration Report By Patrick Gloster Page 1 Patrick Gloster Contents Section 1 – Accessing code and results………………………………………....Page 4 Section 2 – Introduction…………………………………………………………Page 5 Section 3 – The theory of least squares minimization…………………………..Page 5 Section 4 – Problem 1 - Finding the position of a target point using measurements from 4 known observer points in 2D……………………………………….……Page 9 Section 5 – Problem 2 - Finding the position of a set of target points and the observer points themselves using a set of measurements in 2D………………………….Page 12 Section 6 – Problem 3 – The bow-tie problem in 2D…………………………..Page 17 Section 7 – Problem 4 – The bow-tie problem in 3D…………………………..Page 28 Section 8 – Saving data………………………………………………………...Page 41 Section 9 – Current Matlab code…………………………………….…………Page 42 Page 2 Patrick Gloster 1)Accessing code and results All of my files have been transferred from my H drive to the computer PPXP435. They are all contained in a folder entitled ‘CSM Summer 2008’ on the C drive. While writing this document I wrote instructions on how to access the results of different programs. Since then I have reorganised the saved files and code. The following should enable all data to be accessed: When Matlab starts up, choose the current directory (found centrally at the top of the screen) to be ‘C:\CSM Summer 2008\MATLAB\Code’ (Note – I’m not sure if you need the C at the beginning) To run a program type the name of the program followed by the input arguments in brackets, separated by commas. Each program is described in detail in this document, and each description tells you exactly how to start running the program, but they all have the same general format. For example, to run the program bowtie5.m with 100 plate 2 positions and save it to a file named ‘bowtie5_10_figure1.bmp’ type ‘bowtie5(100,10)’ in the Matlab command window and press return To look at the code click on the open folder icon, and then choose the Code folder, followed by whichever program you wish to see the code for. All programs with automatic saves will save files to the Matlab Code folder. If you haven’t moved the saved files from this folder after you generate results, they can be accessed as written in the text Previous results have been moved into separate folders within Matlab. To access the text files, first switch directory to ‘C:\CSM Summer 2008\MATLAB’, then type ‘type FOLDER/FILENAME.TXT’ into the command window. For example, the results of bowtie 5 are in the bowtie5 folder. To access the information from the 102nd run, make sure you are in the Matlab directory, then type ‘type bowtie5/bowtie5_102.txt’ Picture files are saved as bmp files and can be accessed by clicking the open folder icon, followed by the folder which has the name of the program they were created by. If the folder looks empty it may be because files of type is set to ‘Matlab files’, in which case change it to ‘all files’ The available folders containing text files are bowtie, bowtie2, bowtie3, bowtie4, bowtie5, bowtie6, bowtie7, bowtie8, bowtiehistogram, newbowtie3, newbowtie4 and newbowtiehistogram. Their names indicate what programs produced the results they contain The other folders contain figures. ‘2D bowtie varying angle’, 2D bowtie varying initial error’, 2D bowtie varying plate range’ and ‘2D bowtie varying platesize’ also contain workspace variables which can be imported into Matlab. These variables are vectors and a matrix containing the results of tests described towards the end of section 5 in this document Currently the results of the 3D programs are NOT saved automatically. I have been saving the histograms manually in designated folders. If you wish to save information automatically, follow the format of lines 204-255 in bowtie5. Also see the section on saving data in this document The folder ‘CSM Summer 2008’ also contains a folder called ‘Excel data’ which has a spreadsheet containing data generated when trying to find the optimum position of plate 2 in the 2D bowtie set-up. This is referenced on page 24 of this write up. Page 3 Patrick Gloster Finally there is a folder named ‘Write ups’, which contains various documents written as the project progressed including odd things I sent to Patrick Brockill, a copy of this write up, and some files which I used to write down some things as they happened. All useful information has been included in this write up, but I have left these other documents intact. Pages 12 and 13 of the file ‘locating observer coordinates’ contain some information on allocating weights to different errors and may come in useful in the future. If the Matlab help section is not helping I suggest using the website: http://www.mathworks.fr/matlabcentral/ Page 4 Patrick Gloster 2. Introduction The goal of this project was to create a program to calculate the positions of one endplate of the compact straightness monitor with respect to the other. One endplate has 4 launch pairs, arranged in a cross. The other has 4 retro-reflectors, also arranged in a cross. The positions are calculated using the method of least squares. The program also had to include a simulation of positions, so that we could compare calculated values with true values of the positions. The project was approached by writing several programs for different versions of the problem, with gradually increasing complexity. For completeness I include a full account of all of the code written, but many programs are no longer used at all. For each section I will point out the programs which are the end products of that section. The reason for so many discarded programs is the diverse approach taken for each version of the problem. For the first section I worked almost exclusively with my own least squares minimization function; however as the project progressed I experimented with various minimization functions already present in Matlab, and finally settled on a non-linear least squares minimization function called lsqnonlin. 3. The Theory of Least Squares Minimization If we have a set of measurements in a vector, L , which fit a set of formulas involving the position variables in the vector X , we can write our equation as A X L where A is some matrix (the design matrix), operating on the position vector X to give us our measurements. Due to the limits on the precision of measuring instruments, taking redundant physical measurements gives inconsistent results. This means no unique solution exists; all we can do is make a unique estimate of the solution. The chosen method is the method of least squares. We can specify an initial approximation to the solution vector, X 0 , which we perturb slightly to give the total solution vector X X 0 X Since all measurements are fundamentally inaccurate, however, the measured lengths L will give inconsistent results. We therefore add a ‘residual’ vector V which will cancel out the inconsistencies. Thus A X L is satisfied where L L V The method of least squares tells us that the best estimate for X is the estimate which minimises the sum of the squares of the inconsistencies, i.e. minimises V T V Page 5 Patrick Gloster The weight matrix Our measurements (which make up the elements of L) will not necessarily all have the same precision. This is accounted for by assigning a known ‘weight’ matrix to each measurement. P is the matrix whose elements are these weights. In this case the least squares criterion becomes V T PV = minimum Linearization of the problem Mathematical models expressing the relationship between the total observed and solution vectors X and L have the general form: F ( X , L) 0 Linearization of these models is accomplished by replacing the nonlinear functions F by their Taylor’s series linear approximation, expanded about the point defined by the initial approximation to the solution vector (X0), and the measured values of the observation vector (L). Remembering that X X 0 X and L L V linearization gives: F ( X X 0 , L L) F ( X X 0 , L L) F ( X , L) F ( X 0 , L) X V 0 X L Or W AX BV 0 F ( X X 0 , L L) F ( X X 0 , L L) Where W F ( X 0 , L) , A and B X L In our problem the observed quantities can be explicitly expressed as functions of the parameters X (i.e. we can write our problem as F ( X ) L ) In this case, linearization gives us F ( X X 0 ) F(X ) L F(X 0 ) X (L V ) 0 X Or W AX V 0 L Where W F ( X 0 ) L Page 6 Patrick Gloster F1 F1 X 1 X 2 F ( X X 0 ) F2 A = X X 1 Our problem is now: V T PV = minimum X’s subject to constraint: W AX V 0 To implement this constraint we use Lagrange multipliers. Lagrange Multipliers Suppose we want to minimize some function f1 ( x, y) subject to some constraint f 2 ( x, y) 0 The Lagrange method consists of 3 steps: 1) Form the variation function f1 ( x, y) kf 2 ( x, y) where k is an undetermined constant called the Lagrange multiplier 2) Set the derivatives of the variation function to zero 0, 0 x y 3) Solve the three equations f 2 ( x, y) 0, 0, 0 x y In our case we have V T PV = minimum X’s subject to constraint: W AX V 0 So we form the variation function Page 7 Patrick Gloster V T PV 2 K T ( AX V W ) Here K is a column vector of Lagrange multipliers, and has been multiplied by a factor of 2 for convenience. The derivatives are: 2V T P 2K T 0 V which becomes PV K 0 and 2K T A 0 X which becomes AT K 0 So the three equations we want to solve are W AX V 0 PV K 0 AT K 0 These equations can be combined into a single hypermatrix equation (a hypermatrix being a matrix whose elements are themselves matrices) P BT 0 V 0 I 0 A K W 0 0 AT 0 X 0 This can be partitioned (see pages 114-116 of Wells and Krakiwsky) to give X ( AT PA) 1 AT PW K P( AX W ) V P 1 K AX W The solution we are interested in here is X. Thus, if we can determine the matrices A, P and W we can solve for X. Page 8 Patrick Gloster 4. Problem 1 – Finding the position of a target point using measurements from 4 known observer points in 2D The first problem was set up as follows: Figure 1 – 4 observer points with measured distances to the target point (x,y) Here we are measuring the distance of a particular point from four ‘observer points’. We know the positions of these four observers. Our goal was to write a simulation which could generate a large number of target points and, for each of these points, generate measurements from the target point to the 4 observer coordinates, which would be inexact (as all physical measurements are). The method of least squares was then used to calculate an approximate position for the target point and the results were represented graphically. I only created one simulation for this problem, and used my own minimization function. I used the weight matrix A= ( x0 x1 ) ( y 0 y1 ) ( x0 x1 ) 2 ( y 0 y1 ) 2 ( x0 x1 ) 2 ( y 0 y1 ) 2 ( x0 x 2 ) ( y0 y2 ) ( x0 x 2 ) 2 ( y 0 y2 ) 2 ( x0 x 2 ) 2 ( y 0 y 2 ) 2 ( x 0 x3 ) ( y0 y3 ) (x x )2 ( y 0 y3 ) 2 ( x 0 x3 ) 2 ( y 0 y 3 ) 2 0 3 ( x0 x1 ) ( y 0 y1 ) ( x0 x 4 ) 2 ( y 0 y4 ) 2 ( x0 x 4 ) 2 ( y 0 y 4 ) 2 x0 and y0 were approximations to the target position, generated in the simulation. You can choose how many points you want to generate and calculate positions of. The program calls several functions while it is running, so what follows is a description of each of these functions, followed by a description of the simulation itself. Page 9 Patrick Gloster F.m The F function takes two vectors describing two positions, and calculates the distance between the two positions. The input must be in vector form. For example, to calculate the distance between the points (2,4) and (0,0) you need to define two vectors X1=[2;4] and X2=[0;0]. By typing ‘F(X1,X2)’ in the command window the distance between these two points will be returned. Locate.m The locate.m function is used to find the location of an unknown point using two measurements and three known locations. It has an input in the Matlab command window of ‘locate(x1,y1,L1,x2,y2,L2,x3,y3)’. These are 3 sets of coordinates describing the 3 known points, and measurements from two of these points to the unknown point. Two possible positions for the target (unknown) point are calculated using the cosine rule. Whichever of these two calculated points is closer to the third known point is selected as the approximate location of the unknown point. As an example, if you have 3 known points, (0,0), (5,3) and (6,8), and you want to know the approximate location of a point which is 2 units from (0,0) and 5 units from (5,3) you would type locate(0,0,2,5,3,5,6,8). One early problem with this function was what happened if there was no real solution. This happened if both of the measurements were less than the true distance. In this case the program returns the point halfway between the two known points as the solution. Blackcircle.m, bluecircle.m, greencircle.m, redcircle.m These functions plot a circle of radius r centred at (x,y) in the colour included in the name. So to plot a red circle centred on (4,2) with a radius of 1, type ‘redcircle(4,2,1) in the Matlab command window. Simulation.m To begin the simulation, type ‘simulation(N)’ in the Matlab command window, where N is the number of simulations you want to run (that is, how many target positions you want to generate and calculate). This starts the program. First a series of GUIs are generated prompting the user to choose the positions of their known observer coordinates. The user is then given the option to jump to a simulation number of their choice, and asked to input the errors on each measurement. This is so that the weight matrix P can be created. There follows a small loop in the code generating a series of random numbers. This is done because the seed of the random number generator used in this program is set to 0. If you want to jump to a particular simulation, the correct number of random numbers needs to be generated beforehand, so that the correct target point is generated. This loop takes care of this by generating the same number of numbers as would be done if you had run through all of the previous simulations rather than jumping to that one. For each simulation the program generates a target position somewhere in a 10 x 10 square. The true distances from each observer point to the target are then calculated, Page 10 Patrick Gloster and measurements are generated by adding a different random number (normally distributed about 0 with a standard deviation of 1) multiplied by the error for that measurement (entered by the user earlier in the program) to each of the measurements. The simulation then generates an approximate position to start the minimization from. This is done using the locate function described above. A figure is plotted showing circles centred on the four observer coordinates, with radii corresponding to the measurements. The minimization proceeds by applying the general formula above; X ( AT PA) 1 AT PW . This is done several times, converging to a solution. At the beginning of the program I define a variable moddx. The iterations of the minimization continue until the modulus of the perturbation to the solution vector drops below moddx. There is a limit on iterations of 1000000 to prevent the program crashing if it is not converging. The user is given the option to plot each successive solution for the simulation, or to jump to the final solution. Figure 2 – Results of the program simulation.m. The user can zoom in on successive iterations Page 11 Patrick Gloster Gauss.m, gauss2.m There were also two functions written to plot histograms, gauss.m and gauss2.m. The first, gauss.m, is identical to the simulation but at the end it calculates the difference between the true and calculated point for every simulation, and then plots a histogram of the differences. To use this function type ‘gauss(N,simulationnumber)’ into the command window, where N is the number of target points you want to generate and simulationnumber is the simulation you want to jump to. The second is more useful, and is run by typing ‘gauss2(N)’ into the command window, where N is again the number of target points you want to generate. The program runs the simulation and then calculates the distance from each observer point to the calculated point. It then calculates the difference between these 4 distances and the 4 true distances. This is done for each simulated target point. After all of the simulations 4 histograms are plotted of the differences between true and calculated distances to each observer point. Figure 3 – An example of 4 histograms of the differences between calculated and true distances, produced by gauss2.m. Notice the bottom two histograms having a wider spread than the top two. This is interesting and something which I could never adequately explain. However, it was felt that rather than dwelling on this problem we should move on to a new set-up. 5. Problem 2 – Finding the position of a set of target points and the observer points themselves using a set of measurements in 2D The next stage in the project was to introduce uncertainty on where the observer points themselves were, and use many measurements to different target points to calculate where these observer points were (and where the target points were). In problem 1, we solved for each simulated target point separately, and had just 2 unknowns (the x and y coordinates of the target point). In this problem, however, we had to solve for all of the simulated target points simultaneously, since it was only by considering all of the sets of measurements at once that we could solve for the positions of the observer points. For example, we could generate 300 target positions; this would give us 600 unknowns (300 x coordinates, 300 y coordinates). We also have the 8 unknown coordinates describing the 4 observer points. Since we can arbitrarily define an origin Page 12 Patrick Gloster and a direction in this problem, we can reduce the number of unknown coordinates describing the 4 observer points to 5 by setting one of the points to be the origin, (0,0), and one to be along the y-axis, so at a position (0,y). In the example we are then left with 605 unknowns. However, for each generated target position we have 4 measurements, one from each of the observer coordinates. This therefore gives us 1200 measurements in total in our example. Since we have more measurements than unknowns we can again use the least squares method to give us a solution for all of the variables. (0,y2) (x3,y3) Figure 4 – The set- up for problem 2. Here we have 3 of the 8 variables set to zero so we only have to work with 5 unknowns (x4,y4) (0,0) At this stage in the project I began to attack the problem on several fronts. Where before I had just used my own minimization function, I now started to use some of Matlab’s inbuilt functions and compare results. This led to many different functions being written, which are described below; while all are of some interest, the primary fucntions which I would use to find solutions to this problem are basicpat.m and its associated objective function somefun.m. Basicpat.m To use this function type ‘basicpat(N)’ in the Matlab command window, where N is the number of target points you wish to generate. This function generates a set of 5 true coordinates describing the positions of the observer points. We only generate 5 because we have a built-in assumption that one of the observer points is at (0,0) and the second is along the y-axis relative to this first point. The function then generates 2N coordinates to describe the N target points. By adding random errors to all of the target and observer coordinates the initial values are obtained to go into the vector X0. The function then calculates the actual distances from each observer point to each target point, and generates our measurements by adding an error to these distances. The error on the measurements (1E-5) is much less than the error used to give us our initial values in X0 (0.3). An anonymous function is then defined which calls a function called somefun, giving it the number of target points, the measurements, and the current values in the vector Page 13 Patrick Gloster X. This anonymous function is necessary because all of Matlab’s inbuilt functions can only have one input, so we define an anonymous function which only has one input (in this case X) and calls another function which has several inputs (in this case X, measurements and N). We then use Matlab’s function to minimize the anonymous function. The function then uses a built-in function called lsqnonlin to minimize the anonymous function and give us a solution for the vector X. We use lsqnonlin since it is a non- linear least squares minimizer. Somefun.m The function called in basicpat.m is somefun.m. This takes the vector X, which contains our current answers for all variables, and works out what measurements we would get if these were the values. It then works out the difference between these and the actual measurements we have. By minimizing this we are therefore minimizing the difference between the true measurements and these ‘predicted’ measurements. Note – lsqnonlin sums and squares the residuals for us, which is why all somefun does is calculate a vector of the differences. Findobservertwod.m Here we come to a program using my own least squares minimization. It can be run by typing ‘findobservertwod(N)’ in the Matlab command window. It solves the same problem as basicpat, with 5 unknown observer coordinates. This was encouraging because it gave the same answers as basicpat, indicating that they were both working. However, it was much slower, which is why we decided to use Matlab’s least squares minimizer instead. In writing this function I needed to write another function which could calculate the symbolic form of the design matrix A. This was done in the function makejacobian.m. Below is a comparison of the results using basicpat and findobservertwod. Note the time taken for each function. Makejacobian.m This function is called in findobservertwod.m, and returns the symbolic form of the matrix A. The values of the vector X are then substituted into the matrix. This program proved difficult to write because I needed to be able to define an arbitrary number of symbolic variables, which is a non-trivial problem in Matlab. Eventually this was solved by asking for help on the internet forum Matlab central. Y2 X3 Y3 X4 Y4 Time taken (seconds) True value 0.219 0.047 0.678 0.679 0.934 - Lsqnonlin 0..2204 0.0462 0.6827 0.6764 0.9406 0.186 calculated value Findobservertwod 0.2204 0.0463 0.6827 0.6764 0.9405 112.5 calculated value Page 14 Patrick Gloster Basicpat2.m This function is identical to the previous function, basicpat, except it uses a different inbuilt minimizer called fminunc. This is not a least squares minimizer though; instead it is equipped with several types of search and therefore, while it is more powerful, it is not necessary to use it in our problem. I experimented with using this function only as a step to using fmincon, a minimization function which can handle constraints. To run this function type basicpat2(N) in the command window where N is the number of target points you want to generate. Somefun2.m Basicpat2 calls a different objective function to basicpat because, unlike lsqnonlin, fminunc does not accept an objective function which returns a vector. Somefun2 therefore returns a scalar which is the sum of the squares of the difference between the ‘predicted’ measurements and the actual measurements. It is this that is minimized. Pat.m This function was an extension of the function basicpat.m. To use it type pat(N,M) in the Matlab command window, where N is the number of targets and M is the number of simulations you want to perform. The program performs the minimization using lsqnonlin, and repeats the entire thing M times. For each simulation it calculates the difference between the true position and the calculated position (e.g. true x3-calculated x3), and also the residuals of the objective function for each coordinate. Histograms are then plotted for all 10 of these, showing the data from all of the simulations. Basicpat3.m I wrote this function as an experiment in using constraints. Instead of including just 5 variables for the coordinates of the observer points, I included 6, and applied the constraint that one of them was 0 as part of the minimization. This was done using a Matlab function called fmincon. The function works in a very similar way to basicpat.m and basicpat2.m, except where they have x2=0 built in as a condition, and therefore don’t work with it as a variable, basicpat3 includes it as a variable but then applies the constraint that x2=0 during the minimization. (x2,y2) (x3,y3) Figure 5 – The set-up for the program basicpat3. We include x2 as an unknown variable and constrain it to be zero during the minimization Page 15 (x4,y4) (0,0) Patrick Gloster This program was a success; provided the true position of x2 was indeed close to zero the answers were correct. Even if the true position of x2 was not zero, the constraint was still applied; however, this led to huge inaccuracies in all of the other results. Since basicpat3 works with 6 variables for the observer points instead of 5, I had to use a different objective function again, so I wrote somefun3.m. I also included a matrix in the program which defines the constraints. This has to be passed to fmincon as well as the function to be minimized. Here we see how the constraints work if we are applying a sensible one (i.e. saying that y2=0 when it is close to 0) but if y2 is nowhere near 0 applying the constraint destroys all accuracy in the calculated values: Y2 X3 Y3 X4 Y4 True 0.0010 0.0470 0.6789 0.9347 0.3835 position Calculated 0 0.0476 0.6717 0.9303 0.3923 position Y2 X3 Y3 X4 Y4 True 0.2190 0.0470 0.6789 0.9347 0.3835 position Calculated 0 0.2238 -0.5201 -0.1776 0.9441 position Somefun3.m This is the objective function for basicpat3. It is different to somefun1 and somefun2 because when working out the predicted measurements it includes the variable x2, where before this was set to zero. So where before the distance from the target point (x,y) to point 2 (x2,y2) was ( x 2 ( y y 2) 2 , somefun3 calculates the distance using ( x x 2) 2 ( y y 2) 2 . Newpat.m This function was written to test the use of fmincon without constraints. It is identical to basicpat2.m, working with just 5 variables for the observer coordinates, only it uses fmincon instead of fminunc. Findobservertwod2.m & findobservertwod3.m Both of these functions are solving the same problem as basicpat3 (so 6 unknown observer coordinates) using my own minimization. Again, they call functions to give them symbolic design matrices (named makejacobian2.m and makejacobian3.m respectively). They differ from findobservertwod.m because they apply constraints. These programs were both successful in calculating answers and applying constraints, but were abandoned due to the length of time they took. Page 16 Patrick Gloster 6. Problem 3 – The bow-tie problem in 2D The bow-tie set up is a 2D version of the real experimental apparatus. We have two launch pairs on a fixed plate, and use them to measure the distances to two retro- reflectors on a second plate. It is called the bow-tie set up because the 4 measurements form a bow-tie shape. By moving around the right hand plate (plate 2) to different positions and obtaining a large number of measurements we can determine the coordinates specifying the locations of both plates. 1m (a2,b2) L2 y 10 cm (x2,y2) L4 x L3 z (a1,b1) L1 (x1,y1) Figure 6 – The bow-tie set up in 2D with point source launch pairs In the diagram above the left hand plate (which contains the launch heads for the lasers) is fixed; the right hand plate (which contains the retroreflectors) can move to the left and right, up and down, and rotate about the z axis. There is an error on the laser measurements of 1 nm. The error on our initial guesses for positions is 1mm. For the first program, I modelled the launch pairs as point sources, which is why we only have two points on the left hand plate. By defining the bottom point as the origin, and the top point as lying along the y-axis, we obtain x1=0, y1=0, x2=0. Distance scales are marked on the diagram. The distance between the two reflector points on the right hand plate is constrained to be the same for every position of this plate. An initial attempt to write a function to solve this problem using my own minimization method was quickly abandoned due to incorrect results. Similarly, when I wrote a function to solve this using fmincon, Matlab’s minimization function which allows the application of constraints, we found that the results were only accurate to approximately one tenth of a millimetre. We also found that increasing the number of target points did not improve the accuracy. We thus decided to reformulate the problem in a way which did not require constraints, so that we could use lsqnonlin. Page 17 Patrick Gloster Reformulation of the problem To reformulate the problem I changed the way I defined the position of the second plate. Where before it was defined using the 4 coordinates x1,y1,x2,y2, it is now defined using midx, midy and alpha, where midx and midy are the coordinates of the centre of the plate, and alpha is the angle of orientation of the plate. Old formulation New formulation (x2,y2) d=Distance between reflectors α (midx,midy) (x1,y1) Figure 7 – The old and new formulations of the bow tie problem Before I worked with 4 unknowns for every plate position, x1, y1, x2, y2, and constrained the distance between these two points to be the same for every position of the plate. With the new formulation, there are 3 unknowns for every plate position, midx, midy and alpha. There is also d, which is one variable for all positions (i.e. there is not a new d for every plate position). From midx, midy, alpha and d we can work out the positions x1, y1, x2 and y2. We can therefore solve for midx, midy and alpha and never have to apply any constraints during the minimization; this is already part of the new formulation since in working out x1, y1, x2 and y2 (which we need to do so that we can calculate ‘predicted’ measurements) we always use the same d. As before there are several similar programs for this problem. Each one adds a little complexity to the set-up, so the final programs are bowtie5.m and new_bowtie5_histogram.m. Bowtie2.m This program can be run by typing ‘bowtie2(N,filenumber)’ in the command window, where N is the number of different plate positions you want to generate and filenumber is the number of the file you want to save the results as. The program has b2 set to be 0.1, and the distance between reflector points to be 0.1. The other coordinates of plate 1 are already set to be 0, so we only have one coordinate to determine on this plate. The program then generates our initial values Page 18 Patrick Gloster for these variables by adding a normally distributed random error to them. The error on these initial values is 1mm. Using random numbers I then generate the midx and midy values of the second plate N times. I also generate an angle of orientation of plate 2 for each position. Using these I calculate the true positions of the reflectors for each plate position. Again, by adding a random error in mm the program generates our initial values for the midpoint and angle of the plate. The program then uses the true values for the positions of the reflectors to calculate the distances to the observer points, and adds a random error to give us the measurements. The error on the measurements is 1nm. It is again necessary to use an anonymous function in this program. The anonymous function calls muchfun2.m, the objective function, which calculates predicted measurements based on the current values of all of the variables. The anonymous function is then minimized using the Matlab minimization function lsqnonlin. To achieve high accuracy on the results I set one of the optimization parameters, TolFun to 10-18. The lsqnonlin function then calculates values for all of the variables – b2 and then midx, midy and alpha for N different plate positions. The program then calculates the difference between the true values and the calculated values for midx, midy and alpha and calculates the standard deviation of each. It also plots a histogram of these differences for each of the three parameters. Finally the program saves all information in a text file called myfile_filenumber.txt, where filenumber is the second number you entered in the brackets in the command window. To access this information you need to type ‘type myfile0filenumber.txt’ in the command window. For example, if you enter bowtie2(50,4) in the command window, to see all of the information saved you need to enter ‘type myfile04.txt’ into the command window. The information saved includes a table of initial values, calculated values and true values for midx, midy and alpha, the exit flag, the number of iterations, the number of plate positions, the real value of b2, the calculated value of b2, the real separation of reflectors and the calculated separation of reflectors. The histograms are also saved under 'myfile0filenumber_figure1.bmp'. These pictures can then be found in the bowtie2 folder within the Matlab folder in ‘CSM Summer 2008’. One problem experienced with saving the pictures is that it seems to save them by taking a print screen. This means that if you are looking at something else while the program runs, and it gets to the point in the program when it should save the histograms, it may end up saving whatever is on the screen at the time. Muchfun2 This is the objective function minimized in bowtie2. It uses the current values of midx, midy and alpha to calculate current values of x1, y1, x2 and y2, and then uses these and the current value of b2 to calculate predicted measurements. Finally it finds the difference between these predicted measurements and the true measurements. This is what is minimized. Page 19 Patrick Gloster Switching to finite sized launch pairs on plate 1 In bowtie2.m the lasers are launched from point sources. The next step was to separate them so that each of the four beams was emitted from a different source, and to allow them to protrude slightly from the surface of the plate. We use a separation of the launch heads of 2cm and a protrusion of 2mm. Notice that in this set-up the cross measurements (L3 and L4) are performed by the inside launch heads. This is something that was changed for the next function. 1m (0,b2) L2 y 10 cm 2cm (x2,y2) L4 (a4,b4) x (a3,b3) L3 2mm z (0,0) L1 (x1,y1) Figure 8 – The bow-tie set up in 2D with separated launch heads Bowtie4.m This program runs in the same way as bowtie2.m, but it includes a3, b3, a4 and b4. Due to the inclusion of this it also needs a different objective function, the aptly named muchfun4.m. As before, to run the program type bowtie4(N,filenumber) in the command window. Information and results are saved as bowtie4_filenumber.txt, and can be accessed by typing ‘type bowtie4_filenumber.txt’ in the Matlab command window. Histograms are saved as 'bowtie4_filenumber_figure1.bmp' and can be found in the Matlab folder in ‘CSM Summer 2008’. new_bowtie4_histogram.m I wrote new_bowtie4_histogram.m so that we could examine the results for the positions of the launch heads more accurately. This is run by typing ‘new_bowtie4_histogram(N,n,filenumber)’. This works by generating the true values for N plate positions, as before, and then generating n sets of initial positions and Page 20 Patrick Gloster measurements. The program then solves for each of these n sets of measurements. It then calculates the differences between true and calculated values for b2, a3, b3, a4, b4, midx, midy and alpha. We have n values for the observer points (a’s and b’s), and n x N values for plate 2 parameters (midx’s, midy’s and alphas). Histograms are then plotted for all of these variables. The histograms are automatically saved as 'new_bowtie4_run(filenumber)_histogram.bmp'. Measuring cross distances with the outside launch heads In the previous programs I measured the cross distances L3 and L4 with the inside launch heads. The next step was to switch to measuring these with the outside launch heads and compare results. 1m (0,b2) L2 y 10 cm 2cm (x2,y2) L4 (a4,b4) x (a3,b3) L3 2mm z L1 (0,0) (x1,y1) Figure 9 – The bow-tie set up in 2D with crossed beams Bowtie5.m This program solves the system with the new set-up. Files are saved as bowtie5_filenumber.txt and histograms are saved as bowtie5_filenumber_figure1.bmp. To compare the results of this with the results of bowtie4, which is exactly the same but with the outer heads doing the straight measurements and the inner heads doing the cross measurements, I set all true values to be identical, used the same seed for the random number generator for both, used the same number of targets for each, and put exactly the same true error on each measurement. I obtained the following results: Page 21 Patrick Gloster Program File Seed Midx standard Midy standard number deviation deviation Bowtie4 72 0 1.94E-6 2.49E-6 Bowtie5 18 0 8.03E-10 7.55E-10 Bowtie4 73 12 8.76E-7 1.15E-6 Bowtie5 19 12 3.06E-9 5.32E-9 Bowtie4 74 678 3.82E-6 5.20E-6 Bowtie5 20 678 2.90E-9 5.52E-9 Bowtie4 75 1008 9.01E-10 2.43E-9 Bowtie5 21 1008 3.78E-9 5.93E-9 Bowtie4 76 66.6 1.94E-10 2.87E-9 Bowtie5 22 66.6 3.28E-9 6.26E-9 We see that while the results for bowtie4 fluctuate wildly, the results for bowtie5 are consistently extremely good. We also notice that the results are much more accurate in this unrealistic case, when the error on each measurement is exactly the same (by which I mean the actual difference between the measurement and the true value, not just the size of the possible error). This was further examined and produced the following results: Errors File number Seed Midx Midy standard Standard deviation deviation Identical 18 0 8.03E-10 7.55E-10 Non-identical 23 0 3.46E-9 4.18E-8 Identical 19 12 3.06E-9 5.32E-9 Non-identical 24 12 6.02E-8 1.23E-7 Identical 20 678 2.90E-9 5.52E-9 Non-identical 25 678 2.10E-8 4.76E-8 Identical 21 1008 3.78E-9 5.93E-9 Non-identical 26 1008 2.30E-8 1.05E-7 Identical 22 66.6 3.28E-9 6.26E-9 Non-identical 27 66.6 1.68E-8 2.12E-8 Page 22 Patrick Gloster Figure 10 - An example of the histograms plotted in bowtie5 The mean and the standard deviation Upon closer inspection it was found that the results given by bowtie5 were not what we would expect; the standard deviation was of the order of nanometres, but the average size of the difference between the true and calculated values was significantly larger than this standard deviation. We found gaussian distributions which were not centred on 0. The table below shows some results; seed is the number given to the random number generator to act as the seed. The error on initial measurements was 1mm. 200 targets were used. File number Seed Mean x Standard deviation on x 73 2 1.87E-9 1.94E-9 74 8 1.52E-7 2.75E-8 75 42 6.10E-7 8.38E-8 76 4096 -5.37E-9 3.40E-9 77 11 5.58E-9 1.87E-9 Page 23 Patrick Gloster It was hypothesised that this was due to poor calculated values for the launch points on plate 1, i.e. b2, a3, b3, a4 and b4. This was supported when we ran the program without any error on these parameters, and without including them as variables to be minimized. When this was done the gaussians were distributed about a mean much closer to zero than the standard deviation: Changing the error on initial values Figure 11 – When we had no uncertainty on the positions of the launch heads on plate 1, we got the distribution we expected For the majority of the project I worked with an error of 1mm on our ‘initial estimates’; that is, the error added to the true positions to give us our starting point (which we fed into the minimization program lsqnonlin) was of the order of 1E-3. When I experimented with this value, it produced some shocking results which I have not yet found an explanation for. In each case the same seed, 0, was used. Again there were 200 targets. File Error on Mean x Standard Mean y Standard number initial deviation deviation values on x on y 113 1E-3 1.28E-10 1.36E-10 7.95E-8 1.52E-8 114 1E-5 -9.93E-9 5.89E-9 -7.30E-7 1.89E-7 115 1E-7 -3.93E-10 1.20E-9 3.54E-8 1.33E-8 Page 24 Patrick Gloster This was done with the expectation of nothing happening, since no matter what position you start from I would have expected the final convergence to be the same. However, previous tests (page 80 of the logbook), which had some implicit errors and are thus not included here, support this trend. Further tests were carried out using 50 repeats (rather than just looking at one sample with the same seed each time), and again showed a large and unexpected change in accuracy. Graphs of these results are saved in the Matlab folders ‘bowtie5’. Finding the optimum range of movement of plate 2 The next step was to use the program bowtie5.m and vary the positions of the launch heads, the angle which the second plate could be orientated at, and the distance the second plate could move in the x and y directions, in an effort to determine the optimum arrangement. Initially I proceeded by looking at 5 samples, each with a different seed, and every time I changed a parameter (i.e. the angle or the distance that plate 2 could have) I recorded various values for each of these 5 samples. The values recorded were a) the difference between the calculated value for b2 and the true value b) the difference between the calculated value for a3 and the true value c) the difference between the calculated value for b3 and the true value d) the difference between the calculated value for a4 and the true value e) the difference between the calculated value for b4 and the true value f) the mean of the differences between the calculated and true midx values g) the standard deviation of the differences between the calculated and true midx values h) the mean of the differences between the calculated and true midy values i) the standard deviation of the differences between the calculated and true midy values This information is saved in EXCEL in a file called OPTIMIZATION TABLES. This test was, however, limited in scope; looking at 5 samples wasn’t enough. I therefore used the program new_bowtie5_histogram.m to look at 200 targets with 50 different random seeds. This program is similar to the other new_bowtie_histogram.m files; for each of the 50 repeats the same ‘true’ 200 plate positions are used, but each time a different random error is added to the positions (to give us our initial estimates) and to the distances (to give us our measurements). Using this method I varied the range of angles the plate could have, the range of movement of the plate and the size of the plates (keeping the plates the same size as each other for each case). The results of these tests are saved as variables and can be imported into the Matlab workspace. The variables from the test of varying the range of motion that plate 2 can have in the x and y direction are saved as varying_range.m. The variables from the various angle test are saved as varying_angle.m and the variables of changing the size of the plates are saved as varying_platesize.m.. I also used this method to investigate changing the the error on the initial estimates; these variables are saved under varying_initial_error.m. Page 25 Patrick Gloster In each of these cases the results are saved in a matrix: range_matrix, angle_matrix, platesize_matrix or error_matrix. Each column of the matrix represents a different parameter (each column represents the same thing in each matrix). They are defined below: Column number Parameter 1 Range (+/- m) 2 Angle (radians) 3 Mean difference between true b2 and calculated value (m) 4 Standard deviation of b2 difference 5 Mean of a3 difference 6 Standard deviation of a3 difference 7 Mean of b3 difference 8 Standard deviation of b3 difference 9 Mean of a4 difference 10 Standard deviation of a4 difference 11 Mean of b4 difference 12 Standard deviation of b4 difference 13 Mean of mean difference between true and calculated x values * 14 Standard deviation of mean x differences 15 Mean of standard deviation of x differences 16 Standard deviation of standard deviation of x differences 17 Mean of mean of y differences 18 Standard deviation of mean of y differences 19 Mean of standard deviation of y differences 20 Standard deviation of standard deviation of y differences 21 Mean of mean of alpha differences 22 Standard deviation of mean of alpha differences 23 Mean of standard deviation of alpha differences 24 Standard deviation of standard deviation of alpha differences 25 Mean of all x differences 26 Standard deviation of all x differences 27 Mean of all y differences 28 Standard deviation of all y differences 29 Mean of all alpha differences 30 Standard deviation of all alpha differences 31 Initial error Page 26 Patrick Gloster * We have these odd values because for each of the 50 repeats, we work out an average difference in both x and y, and we also calculate the standard deviation of these differences. In the matrices we therefore include the averages and standard deviations of all of these averages and standard deviations, and also include the mean and standard deviation of all of the values across all of the repeats. Graphs plotting parameters 3-30 against range, angle and initial error are saved as Matlab figures e.g. to see how the mean difference between the true and calculated values for b2 changed with the amount that the second plate could move see the graph entitled ‘b2_mean_vs_range.fig’. The results show a general trend of increasing accuracy with increasing angle, range and plate size. As we continue to increase these parameters still further the accuracy decreases sharply, but I believe this is due to the optimization no longer being stopped by TolFun; instead it is halted by TolX. It is difficult to make specific statements about which values of these parameters would be the best to use because increasing angle from 10 degrees to 20 degrees could for example make the accuracy in working out a2 better, but worsen the accuraracy of b4. Sparsity pattern When implementing the new_bowtie5_histogram.m program it could be very slow due to the 50 repetitions, and the fact that we were varying the parameters to values which led to extremely long convergence times. This prompted me to add an ammendement to the code which returns the jacobian for the first convergence. All of the non-zero elements of this jacobian are then converted to ones, giving us a matrix representing the sparsity pattern of the jacobian. In every subsequent use of lsqnonlin this sparsity pattern is handed to the function. Using the sparsity pattern makes convergence approximately 10 times faster, thus although the first convergence is slow, the rest of the program runs very quickly. Page 27 Patrick Gloster 7. Problem 4 – The bow-tie problem in 3D Plate 1 is modelled as having 8 launch heads: Looking at plate 1 from plate 2 (a4,b4,c4) (a3,b3,c3) (a6,b6,c6) (a8,b8,c8) (a5,b5,c5) (a7,b7,c7) (a2,b2,c2) (a1,b1,c1) I define a1=0, b1=0, c1=0, a4=0, c4=0, c5=0, giving us a plane in z=0. Plate 2 lies in the positive z direction: 10cm Plate 1 Plate 2 z Page 28 Patrick Gloster Plate 2 is modelled as having 4 retroreflectors: Looking at plate 2 from (x2,y2,z2) = point 2 plate 1 d1 (x4,y4,z4) = point 4 d2 (x3,y3,z3) = point 3 z y x (x1,y1,z1) = point 1 Plate 2 can rotate around 3 axes: 1) Rotation about the x-axis where the angle alpha is the angle between the plate and the y-axis y x y α z z 2) Rotation about the y-axis where beta is the angle between plate 2 and the x-axis x x y β z z 3) Rotation about the z-axis where gamma is the angle between plate 2 and the x-axis y x y γ z x Page 29 Patrick Gloster If we define the midpoint of plate 2 as (midx, midy, midz), then using angles alpha, beta and gamma, and the distances d1 and d2, we can define the locations of the 4 retroreflectors: x1 = midx + (d1/2).cosα.sinγ y1 = midy – (d1/2).cosα.cosγ z1 = midz + (d1/2).sinα x2 = midx – (d1/2).cosα.sinγ y2 = midy + (d1/2).cosα.cosγ z2 = midz – (d1/2).sinα x3 = midx + (d2/2).cosβ.cosγ y3 = midy + (d2/2).cosβ.sinγ z3 = midz – (d2/2).sinβ x4 = midx – (d2/2).cosβ.cosγ y4 = midy – (d2/2).cosβ.sinγ z4 = midz + (d2/2).sinβ We define our 8 measurements as follows: L1 ( x1 a 2) 2 ( y1 b2) 2 ( z1 c 2) 2 L2 ( x1 a 4) 2 ( y1 b4) 2 ( z1 c 4) 2 L3 ( x 2 a3) 2 ( y 2 b3) 2 ( z 2 c3) 2 L4 ( x 2 a1) 2 ( y 2 b1) 2 ( z 2 c1) 2 L5 ( x3 a7) 2 ( y3 b7) 2 ( z3 c7) 2 L6 ( x3 a5) 2 ( y3 b5) 2 ( z3 c5) 2 L7 ( x4 a6) 2 ( y 4 b6) 2 ( z 4 c6) 2 L8 ( x 4 a8) 2 ( y 4 b8) 2 ( z 4 c8) 2 Remembering that a1=0, b1=0, c1=0, a4=0, c4=0 and c5=0. Page 30 Patrick Gloster threeDbowtie.m The program threeDbowtie.m is the program to solve for the positions of the plates in three dimensions. It is very much an extension of the 2D programs; where before I defined 4 points on plate 1 and 2 on plate 2, I now define 8 points on plate 1 and 4 on plate 2. Again the objective function calculates the difference between the predicted measurements and the simulated measurements, and is minimized. It is run by entering ‘threeDbowtie(N,n)’ in the Matlab command window, where N is the number of plate 2 positions and n is the number of times you want to generate these points and minimise them. In previous programs which included repeats the same N true positions were used in all of the repeats; it was only the initial estimates and the measurements which changed. In this program, however, an entirely new set of N plate 2 positions is generated. To look at one run you just set n=1. NOTE – An important point to make is that the variables rangeinmidX, rangeinmidY, rangeinmidZ are the amount that can be added or taken away from the starting (x,y,z) value 1m away from the centre of plate 1. So if rangeinmidZ is set to 0.1m the middle of plate 2 can vary from 0.9m away from the centre of plate 1 to 1.1m away from the centre of plate 1. This is the same in the following programs as well. The manner in which the positions of plate 2 are generated is as follows: 1) We take a point 1m away from the centre of plate 1 along the Z-axis 2) We add a random number between -1 and 1 multiplied by the variables rangeinmidX/Y/Z to this point to give us the centre of plate 2 ThreeDbowtie.m minimizes all of the variables at once, and plots histograms of the differences between the true and calculated values. The next step is to do run the minimization in two steps; first we find values for the positions of the launch heads on plate 1, then we use these calculated values as fixed positions and can calculate the positions of plate 2 with a systematic error. It should be noted that for this program and for calibrationbowtie.m I have an error on initial estimates of angle of 1 degree. This can be changed in line 19 of the m-file. calibrationbowtie3.m This program takes a set of N plate 2 positions and calculates the positions of the launch heads on plate 1. These are then used for the rest of the simulation. A further M plate 2 positions are generated and the program calculates the values of midx, midy, midz, alpha, beta and gamma for these positions, always using the same calculated positions for the a’s, b’s and c’s. This program is run by typing ‘calibrationbowtie3(N,M) in the Matlab command window, where N is the number of plate positions you want to use to determine the calibration constants and M is the number of plate positions you want to calculate the positions of using these calibration constants. My initial set up of the launch heads on plate 1 was completely symmetric. This gave results of: Page 31 Patrick Gloster midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean 1.70E-6 -7.82E-6 1.63E-6 -2.75E-6 -8.34E-6 1.52E-4 Standard 3.59E-6 2.62E-6 5.66E-8 7.22E-5 8.34E-5 1.06E-3 deviation We note that gamma is particularly bad. The positions of the launch heads making up one of the launch pairs was then lowered 5 mm relative to the other pair, so where before b5=0.05, b6=0.05, b7=0.05 and b8=0.05 we now had b5=0.045, b6=0.045, b7=0.05 and b8=0.05. This gave the following results: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean 1.29E-7 -3.49E-7 6.09E-8 -3.67E-7 1.16E-8 8.76E-7 Standard 4.08E-8 4.88E-8 1.10E-9 2.56E-6 2.44E-6 3.22E-6 deviation And when repeated: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean -7.50E-8 -2.46E-7 4.60E08 8.51E-10 4.91E-8 1.53E-6 Standard 3.96E-8 3.79E-8 9.98E-10 1.84E-6 1.77E-6 3.19E-6 deviation This indicated that removing the symmetry greatly improved the accuracy of the calibrations. Calibrations.m The next target was to find out how many targets were required to give us a consistent set of calibrations. This was done using the program calibrations.m, which calculates the calibration constants 50 times and plots histograms of the difference between the true and calculated values.To run calibrations.m type ‘calibrations(N)’ into the Matlab command window where N is the number of plate 2 positions you wish to generate. Histograms are saved in the folder ‘Calibration histograms’ in the Matlab folder. Page 32 Patrick Gloster 100 plate 2 positions: a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference Mean -2.62E-8 -6.80E-9 -9.73E-8 -6.97E-8 4.26E-8 9.81E-8 Standard 8.72E-8 9.08E-8 3.74E-7 2.21E-7 2.42E-7 3.77E-7 deviation b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference Mean 3.92E-8 1.53E-7 1.88E-7 9.14E-8 1.00E-7 8.13E-8 9.18E-8 Standard 1.70E-7 6.03E-7 7.51E-7 3.61E-7 3.79E-7 3.56E-7 3.53E-7 deviation c2_difference c3_difference c6_difference c7_difference c8_difference Mean -1.76E-8 -1.75E-8 -1.80E-8 -1.71E-8 9.31E-10 Standard 6.72E-8 7.30E-8 6.97E-8 7.05E-8 6.01E-9 deviation 200 plate 2 positions: a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference Mean 7.20E-9 1.21E-8 4.67E-9 8.57E-9 7.77E-9 2.25E-9 Standard 6.05E-8 6.06E-8 1.04E-7 7.37E-8 7.72E-8 1.06E-7 deviation b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference Mean 1.46E-8 8.59E-9 4.38E-9 -5.72E-11 1.30E-8 7.02E-9 2.66E-10 Standard 5.98E-8 1.70E-7 2.06E-7 1.01E-7 1.25E-7 1.09E-7 9.69E-8 deviation c2_difference c3_difference c6_difference c7_difference c8_difference Mean 1.02E-9 1.43E-10 8.10E-10 2.53E-10 3.12E-10 Standard 1.83E-8 2.01E-8 1.92E-8 1.94E-8 4.06E-9 deviation Page 33 Patrick Gloster 300 plate 2 positions: a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference Mean -1.15E-9 4.22E-9 -7.71E-9 -3.34E-9 8.13E-9 7.24E-9 Standard 4.72E-8 5.73E-8 5.48E-8 4.12E-8 4.11E-8 5.36E-8 deviation b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference Mean 8.67E-9 1.55E-8 2.24E-8 1.07E-8 1.41E-8 9.88E-9 1.01E-8 Standard 4.22E-8 9.45E-8 1.13E-7 5.50E-8 6.44E-8 8.62E-8 5.40E-8 deviation c2_difference c3_difference c6_difference c7_difference c8_difference Mean -1.12E-9 -1.93E-9 -8.32E-10 -2.27E-9 -9.09E-10 Standard 1.04E-8 1.03E-8 9.94E-9 1.08E-8 2.87E-9 deviation 400 plate 2 positions: a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference Mean -2.02E-9 -4.88E-9 5.54E-10 1.15E-9 -5.73E-9 1.66E-9 Standard 3.49E-8 4.62E-8 4.19E-8 3.71E-8 4.29E-8 4.27E-8 deviation b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference Mean 6.70E-9 -3.24E-9 4.63E-10 -4.34E-10 1.38E-9 3.86E-9 3.42E-11 Standard 3.01E-8 6.78E-8 8.47E-8 3.95E-8 4.69E-8 5.75E-8 3.99E-8 deviation c2_difference c3_difference c6_difference c7_difference c8_difference Mean 1.71E-10 4.57E-10 4.31E-10 1.77E-10 2.42E-10 Standard 7.48E-9 8.68E-9 7.76E-9 8.64E-9 2.66E-9 deviation 500 plate 2 positions: a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference Mean -5.80E-9 -4.99E-10 -4.83E-10 1.88E-9 -4.10E-9 1.58E-9 Standard 3.70E-8 3.76E-8 4.35E-8 3.73E-8 3.04E-8 4.20E-8 deviation b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference Mean 2.53E-9 -4.95E-9 -1.25E-9 -1.32E-9 2.05E-9 -3.99E-9 -1.36E-9 Standard 2.97E-8 6.46E-8 7.83E-8 3.90E-8 5.44E-8 4.97E-8 3.88E-8 deviation c2_difference c3_difference c6_difference c7_difference c8_difference Mean 3.04E-10 1.13E-10 5.84E-11 2.99E-10 7.18E-11 Standard 6.99E-9 7.55E-9 7.36E-9 7.09E-9 2.36E-9 deviation Page 34 Patrick Gloster When attempting to do 600 plate 2 positions Matlab gave an error ‘Out of memory’. We can see that the accuracy of the locations of the launch heads is increasing with an increasing number of plate 2 positions, so it would be advisable to use as many plate 2 positions as possible. However, time is an issue. When using calibrationbowtie3 to calculate the calibration constants using 500 targets, and then using these calculated values to calculate the position of 100 plate 2 positions, the time taken was 416 seconds. Introducing further offsets on the launch heads for calibrationbowtie3.m We saw earlier that lowering one of the launch pairs with respect to its opposite results in a far better accuracy when calculating the positions of plate 2 due to the removal of the symmetry. Subsequent testing should include adding further offsets, with launch heads sticking out of the plate varying amounts and different distances to the centre of the plate for each launch head. The previous settings for calibrationbowtie3.m were: a1 0 b1 0 c1 0 a2 0 b2 0.02 c2 0 a3 0 b3 0.08 c3 0 a4 0 b4 0.1 c4 0 a5 -0.05 b5 0.045 c5 0 a6 -0.03 b6 0.045 c6 0 a7 0.03 b7 0.05 c7 0 a8 0.05 b8 0.05 c8 0 giving results: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean 1.29E-7 -3.49E-7 6.09E-8 -3.67E-7 1.16E-8 8.76E-7 Standard 4.08E-8 4.88E-8 1.10E-9 2.56E-6 2.44E-6 3.22E-6 deviation And when repeated: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean -7.50E-8 -2.46E-7 4.60E08 8.51E-10 4.91E-8 1.53E-6 Standard 3.96E-8 3.79E-8 9.98E-10 1.84E-6 1.77E-6 3.19E-6 deviation Page 35 Patrick Gloster The next test used the following settings for the launch head positions: a1 0 b1 0 c1 0 a2 0 b2 0.02 c2 0 a3 0 b3 0.08 c3 0 a4 0 b4 0.1 c4 0 a5 -0.05 b5 0.045 c5 0 a6 -0.03 b6 0.045 c6 0 a7 0.028 b7 0.05 c7 0 a8 0.048 b8 0.05 c8 0 I have highlighted in red the alterations. We have shifted the two launch pair on the right of the diagram 2mm towards the centre. This means the launch heads on the hoizontal are no longer equidistant from the centre. Looking at plate 1 from plate 2 (a4,b4,c4) (a3,b3,c3) (a6,b6,c6) (a8,b8,c8) (a5,b5,c5) (a7,b7,c7) (a2,b2,c2) (a1,b1,c1) The results for this set up were: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean -1.86E-9 -1.07E-7 1.47E-8 -2.43E-8 2.01E-8 4.22E-7 Standard 1.80E-8 1.48E-8 5.09E-10 6.44E-7 6.91E-7 1.67E-7 deviation Figure saved as ‘figure1’ in the calibrationbowtie3 folder in the Matlab folder. When repeated: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Page 36 Patrick Gloster Mean -7.91E-9 -6.22E-8 1.05E-8 -6.14E-9 2.21E-8 4.26E-7 Standard 1.50E-8 1.56E-8 4.19E-10 5.37E-7 4.27E-7 1.58E-6 deviation Figure saved as ‘figure2’ in the calibrationbowtie3 folder. We see that this appears to have increased accuracy further when compared to the previous results shown on page 34. The next test was to introduce an asymmetry in the vertical components on plate 1. I chose to shift the upper launch pair closer to the centre of the plate. However, I did not try intoducing a horizontal shift because the top launch head (the launch head at (a4, b4, c4)) is defined to be at (0, b4, 0). While it is possible to change this, and worth testing in the future, this must be done carefully because it would involve rewriting the objective function; this is because currently the objective function has the value of a4=0 hardwired into it. The new value would have to replace this in the objective function if we varied a4. The new set-up was therefore: a1 0 b1 0 c1 0 a2 0 b2 0.02 c2 0 a3 0 b3 0.078 c3 0 a4 0 b4 0.098 c4 0 a5 -0.05 b5 0.045 c5 0 a6 -0.03 b6 0.045 c6 0 a7 0.028 b7 0.05 c7 0 a8 0.048 b8 0.05 c8 0 This produced the following results: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean 1.77E-8 -9.52E-8 5.28E-9 -1.25E-7 5.68E-8 4.99E-7 Standard 1.43E-8 1.46E-8 6.42E-10 1.43E-7 1.90E-7 1.28E-6 deviation Figure saved as ‘figure3’ in the calibrationbowtie3 folder in the Matlab folder. When repeated: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean -1.86E-9 -1.07E-7 1.47E-8 -2.43E-8 2.01E-8 4.22E-7 Standard 1.80E-8 1.48E-8 5.09E-10 7.44E-7 6.91E-7 1.67E-6 deviation Page 37 Patrick Gloster Figure saved as ‘figure4’ in the calibrationbowtie3 folder in the Matlab folder. This has not had any obvious effects. It should be noted that these are preliminary tests to give us an idea of any change; rigorous testing would involve the use of calibrationbowtie4 with these difference settings for a’s, b’s and c’s. The next test was introducing some protrusion from the plate: a1 0 b1 0 c1 0 a2 0 b2 0.02 c2 0 a3 0 b3 0.078 c3 0 a4 0 b4 0.098 c4 0 a5 -0.05 b5 0.045 c5 0 a6 -0.03 b6 0.045 c6 0 a7 0.028 b7 0.05 c7 0.002 a8 0.048 b8 0.05 c8 0.002 midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean 3.69E-9 -4.66E-9 -3.66E-9 -2.74E-8 6.35E-8 -2.10E-7 Standard 1.25E-8 1.29E-8 4.80E-10 1.63E-7 1.59E-7 1.98E-6 deviation Figure saved as ‘figure5’ in the calibrationbowtie3 folder in the Matlab folder. When repeated: midX midY midZ Alpha Beta Gamma difference difference difference difference difference difference Mean 2.19E-8 -6.18E-9 4.56E-9 3.96E-8 2.82E-8 5.65E-9 Standard 1.28E-8 1.20E-8 4.49E-10 2.01E-7 2.57E-7 1.02E-6 deviation Figure saved as ‘figure6’ in the calibrationbowtie3 folder in the Matlab folder. Again, there are no huge differences and this needs to be investigated using calibrationbowtie4. Time issues Page 38 Patrick Gloster A problem I am currently experiencing with my programs is the amount of time it takes to run each one. Since more plate 2 positions corresponds to a better calibration, we want to use a large number of them. However, even using just 300 targets takes 131 seconds per calibration. This means looking at 50 calibrations using calibrationbowtie4.m will be incredibly time consuming. If it is possible to speed up the calibration stage this would be extremely beneficial. Due to these problems with the amount of time taken for each calibration, I am looking at calibrationbowtie4 with just 300 plate 2 positions for the calibration, and a further 100 plate 2 positions calculated with that set of calibration constants, for 50 different calibrations. In the future we would like to have more plate 2 positions for the calibration, but this is the largest number of positions which I can do in the time available. The true positions of the launch heads on plate 1 were set up as follows: a1 0 b1 0 c1 0 a2 0 b2 0.02 c2 0 a3 0 b3 0.08 c3 0 a4 0 b4 0.1 c4 0 a5 -0.05 b5 0.045 c5 0 a6 -0.03 b6 0.045 c6 0 a7 0.03 b7 0.05 c7 0 a8 0.05 b8 0.05 c8 0 We obtained the following results: midX midY midZ alpha beta gamma mean of 3.02E-9 -2.80E-8 3.99E-9 -3.10E-9 -1.00E-8 -4.04E-8 mean standard 4.42E-8 6.12E-8 9.42E-9 6.71E-8 8.75E-8 4.37E-7 deviation of mean mean of 1.52E-8 1.56E-8 5.79E-10 3.39E-7 3.32E-7 1.72E-6 standard deviation standard 2.41E-9 2.77E-9 1.60E-10 3.02E-7 2.78E-7 3.30E-7 deviation of standard deviation Figures are saved in the folder calibrationbowtie4 as means1 and SD1. If we look at the mean and standard deviation of the standard deviations we can see that once we have calibrated plate 1 the difference between the calculated and true values has a very consistent standard deviation of less than a nanometre on midZ and Page 39 Patrick Gloster about 15 or 16nm on midX and midY. The average disparity between the true and calculated values is approximately 30nm on midX and midY, and 4nm on midZ. I find this encouraging since only 300 positions were used for the calibration and using 400 or 500 targets has been shown to give better calibrations. I also did a quick check using only 200 targets for each calibration, with the same true positions of the launch heads as above, and obtained the following results: midX midY midZ alpha beta gamma mean of 2.89E-7 -9.30E-7 1.16E-7 -5.01E-7 5.63E-7 1.23E-4 mean standard 2.37E-6 3.90E-6 6.75E-7 3.40E-6 4.22E-6 2.14E-4 deviation of mean mean of 5.69E-7 6.29E-7 1.58E-8 2.00E-5 2.08E-5 1.88E-3 standard deviation standard 8.00E-7 1.20E-6 2.20E-8 2.71E-5 3.05E-5 1.88E-3 deviation of standard deviation Figures are saved as means2 and SD2 in calibrationbowtie4. We can see here the huge difference made by using just 200 plate 2 positions for the calibration instead of 300. Page 40 Patrick Gloster 8) Saving data Any programs I have written which record data automatically save two types: pictures of the histograms and text files containing different variables calculated/present in the function. Saving figures When the figure window is created a figure handle is assigned to it by h1=figure('Name',Whatever name you want’), where the figure handle is h1. This is defined BEFORE the figure is plotted The figure is then plotted by the command hist(vector, number of bins) The figure is saved by writing ‘saveas(figure handle, filename)’ e.g. ‘saveas(h1,myfigure)’ To avoid overwriting saved files many of my functions have an input argument which is a number assigned to the filename you save the picture as. For example in bowtie5 the figures are saved by ‘saveas(h1,sprintf('bowtie5_%d_figure1.bmp',andy));’ where sprintf is used to read the number represented by andy and put this in the place of the ‘%d’. If andy was 10, the histogram would be saved as bowtie5_10_figure1.bmp Warning – this method of automatic saving appears to be a fairly poor method, since it seems to be some form of print screen and sometimes it takes a picture of some other application open at the time, rather than the picture you want saved To access the pictures just open the Matlab folder. If the program you were running is saved in the Code folder (as it probably will be), the figure will be automatically saved in this folder. Double click on it to open the figure. Saving text files First we open a file to save the information to with the line fid = fopen(filename, 'w'); where the w tells Matlab we want are opening the file with the intention of writing information to it We then write information using ‘fprintf(fid,’words = %d’, variable name)’. As an example, to write the word ‘midX’ next to the value for ‘midX’ we would have the line ‘fprintf(fid,’midX = %d’, midX)’ Use \n for a new line To close the file we have the line fclose(fid); To access the text files which are in a particular folder make sure you have selected that folder in the current directory at the top centre of the screen then type ‘type filename.txt’ in the Matlab command window. If you are running programs saved in the Matlab/Code folder then you will already be in this directory and all text files will be automatically saved here. If you haven’t moved them then just typing ‘type filename.txt’ will open them Page 41 Patrick Gloster If the files are in a different folder and you do not want to change directory type ‘type foldername/filename.txt’ to open them 9) Current Matlab code This section contains copies of the 3 most recent programs I have written, which will be used in the future. To change the range of movement of plate 2 alter platerangeinX, platerangeinY and platerangeinZ. All distance units are in metres, and all angle units are in radians. Just to reiterate, the variables rangeinmidX, rangeinmidY and rangeinmidZ define how far in either direction from a central point the plates can move along that axis, so if these are all set to 0.1 then the full range of the midpoints of the plate will be 0.2m3 about a central point. The manner in which the positions of plate 2 are generated is as follows: 3) We take a point 1m away from the centre of plate 1 along the Z-axis 4) We add a random number between -1 and 1 multiplied by the variables rangeinmidX/Y/Z to this point to give us the centre of plate 2 CALIBRATIONBOWTIE3.M function calibrationbowtie3(N,M) format long %This sets the seed for the random number by the clock time=clock; seed=time(6)+time(5)+time(4)+time(3)+time(2)+time(1); rand('seed',seed); randn('seed',seed); %These are the distances between retroreflectors on plate 2 d1=0.1; d2=0.1; initialerror=1E-3; %This is the error on our initial estimates of the positions, 1mm lasererror=1E-9; %This is the error on the laser measurements, 1nm distancebetweenplates=1; %This is the rough distance between the plates along the Z axis, 1 metre angleerror=(pi/180); %This is the error on how well we can determine the angle, currently 1 degree %These are the ranges of movement along each axis that plate 2 has platerangeinX=0.01; platerangeinY=0.01; platerangeinZ=0.01; %These are the sizes of the angles that plate 2 can rotate through alphavariation=2.*pi/9; betavariation=2.*pi/9; gammavariation=2.*pi/9; %Define the points on plate 1 Page 42 Patrick Gloster a1=0; b1=0; c1=0; a2=0; b2=0.02; c2=0; a3=0; b3=0.078; c3=0; a4=0; b4=0.098; c4=0; a5=-0.05; b5=0.045; c5=0; a6=-0.03; b6=0.045; c6=0; a7=0.028; b7=0.05; c7=0.002; a8=0.048; b8=0.05; c8=0.002; %Generate the midpoint values of plate 2 for N positions midx=platerangeinX.*2.*(rand(N,1)-0.5); midy=(b4/2)+platerangeinY.*2.*(rand(N,1)-0.5); midz=distancebetweenplates+platerangeinZ.*2.*(rand(N,1)-0.5); %Generate the angles of orientation of plate 2 for each of N positions alpha=2.*(rand(N,1)-0.5).*alphavariation; beta=2.*(rand(N,1)-0.5).*betavariation; gamma=2.*(rand(N,1)-0.5).*gammavariation; %Generate the true positions of points 1, 2, 3 and 4 x1=midx+(d1/2).*cos(alpha).*sin(gamma); y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); Page 43 Patrick Gloster x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); truevalues=[midx,midy,midz,alpha,beta,gamma]; truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4]; truevalues=truevalues'; truevalues=truevalues(:); truepositions=truepositions'; truepositions=truepositions(:); %Generate our starting approximations for a's, b's and c's. Remember, a1=0, %b1=0, c1=0, a4=0, c4=0 and c6=0 so aren't included as variables approx_a2=a2+randn.*initialerror; approx_b2=b2+randn.*initialerror; approx_c2=c2+randn.*initialerror; approx_a3=a3+randn.*initialerror; approx_b3=b3+randn.*initialerror; approx_c3=c3+randn.*initialerror; approx_b4=b4+randn.*initialerror; approx_a5=a5+randn.*initialerror; approx_b5=b5+randn.*initialerror; approx_a6=a6+randn.*initialerror; approx_b6=b6+randn.*initialerror; approx_c6=c6+randn.*initialerror; approx_a7=a7+randn.*initialerror; approx_b7=b7+randn.*initialerror; approx_c7=c7+randn.*initialerror; approx_a8=a8+randn.*initialerror; approx_b8=b8+randn.*initialerror; approx_c8=c8+randn.*initialerror; %Generate our initial approximations to d1 and d2: approx_d1=d1+randn.*initialerror; approx_d2=d2+randn.*initialerror; %Generate our approximations to the midpoints of plate 2 approx_midx=midx+randn(N,1).*initialerror; approx_midy=midy+randn(N,1).*initialerror; Page 44 Patrick Gloster approx_midz=midz+randn(N,1).*initialerror; %Generate our approximations to the angles of orientation of plate 2 approx_alpha=alpha+(randn(N,1).*angleerror); approx_beta=beta+(rand(N,1).*angleerror); approx_gamma=gamma+(rand(N,1).*angleerror); approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro x_beta,approx_gamma]; approx_values=approx_values'; approx_values=approx_values(:); %Calculate the true distances D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); D4=sqrt((x2).^2+(y2).^2+(z2).^2); D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); %Generate the measurements by adding a random error to each distance L1=D1+randn(N,1).*lasererror; L2=D2+randn(N,1).*lasererror; L3=D3+randn(N,1).*lasererror; L4=D4+randn(N,1).*lasererror; L5=D5+randn(N,1).*lasererror; L6=D6+randn(N,1).*lasererror; L7=D7+randn(N,1).*lasererror; L8=D8+randn(N,1).*lasererror; measurements=[L1,L2,L3,L4,L5,L6,L7,L8]; %Form a vector of our initial estimates of all variables X0=[approx_a2;approx_b2;approx_c2;approx_a3;approx_b3;approx_c3;appro x_b4;... approx_a5;approx_b5;approx_a6;approx_b6;approx_c6;approx_a7;approx_b7 ;... approx_c7;approx_a8;approx_b8;approx_c8;approx_d1;approx_d2;approx_va lues]; %We need to use an anonymous function as the objective function because %superhappyfun has 2 input arguments and lsqnonlin can only minimize %functions with 1 input argument anon = @(X) superhappyfun(measurements,X); %Set lsqnonlin to do one iteration. This will quickly perform the %minimzation with just one iteration, giving us a result with very bad %accuracy. However, what we are interested in is the jacobian which is Page 45 Patrick Gloster %returned. Using this jacobian the sparsity structure is obtained and then %handed to lsqnonlin for the next minimization, which has options suitable %for our minimization and will run much faster using the sparsity %structure. options=optimset('MaxIter',1); [dummy_X,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,... dummy_lambda,dummy_jacobian]=lsqnonlin(anon,X0,[],[],options); %Derive the sparsity structure jac=dummy_jacobian; for a=1:8.*N for b=1:6.*N+20 if jac(a,b)~=0 jac(a,b)=1; end end end %Set the options for lsqnonlin suitable for a good minimization options=optimset('TolFun',1E-18,'TolX',1E- 10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac obPattern',jac); %Use the minimization function lsqnonlin to calculate values for %the positions of the launch heads on plate 1. These values are only %calculated once [X,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(anon,X 0,[],[],options); calibrations=X(1:20); %Clear our previous variables. We can now calculate the positions of a set %of plate 2 positions using the calibration constants we have just %calculated clear midx midy midz alpha beta gamma x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4 clear truevalues truepositions approx_midx approx_midy approx_midz clear approx_alpha approx_beta approx_gamma approx_values clear D1 D2 D3 D4 D5 D6 D7 D8 L1 L2 L3 L4 L5 L6 L7 L8 measurements %Now generate M plate 2 positions and, using the calibrations we have just %worked out, calculate their positions %Generate the midpoint values of plate 2 for M positions midx=platerangeinX.*2.*(rand(M,1)-0.5); midy=(b4/2)+platerangeinY.*2.*(rand(M,1)-0.5); midz=distancebetweenplates+platerangeinZ.*2.*(rand(M,1)-0.5); Page 46 Patrick Gloster %Generate the angles of orientation of plate 2 for each of M positions alpha=2.*(rand(M,1)-0.5).*alphavariation; beta=2.*(rand(M,1)-0.5).*betavariation; gamma=2.*(rand(M,1)-0.5).*gammavariation; %Generate the true positions of points 1, 2, 3 and 4 x1=midx+(d1/2).*cos(alpha).*sin(gamma); y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); truevalues=[midx,midy,midz,alpha,beta,gamma]; truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4]; truevalues=truevalues'; truevalues=truevalues(:); truepositions=truepositions'; truepositions=truepositions(:); %Generate our approximations to the midpoints of plate 2 approx_midx=midx+randn(M,1).*initialerror; approx_midy=midy+randn(M,1).*initialerror; approx_midz=midz+randn(M,1).*initialerror; %Generate our approximations to the angles of orientation of plate 2 approx_alpha=alpha+(randn(M,1).*angleerror); approx_beta=beta+(rand(M,1).*angleerror); approx_gamma=gamma+(rand(M,1).*angleerror); approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro x_beta,approx_gamma]; approx_values=approx_values'; approx_values=approx_values(:); %Calculate the true distances D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); D4=sqrt((x2).^2+(y2).^2+(z2).^2); D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); Page 47 Patrick Gloster D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); %Generate the measurements by adding a random error to each distance L1=D1+randn(M,1).*lasererror; L2=D2+randn(M,1).*lasererror; L3=D3+randn(M,1).*lasererror; L4=D4+randn(M,1).*lasererror; L5=D5+randn(M,1).*lasererror; L6=D6+randn(M,1).*lasererror; L7=D7+randn(M,1).*lasererror; L8=D8+randn(M,1).*lasererror; measurements=[L1,L2,L3,L4,L5,L6,L7,L8]; X00=approx_values; anon2 = @(XX) calibrationfun(measurements,XX,calibrations); %Obtain the sparsity structure for this next minimization options=optimset('MaxIter',1); [dummy_XX,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,.. . dummy_lambda,dummy_jacobian2]=lsqnonlin(anon2,X00,[],[],options); jac2=dummy_jacobian2; for a=1:8.*M for b=1:6.*M if jac2(a,b)~=0 jac2(a,b)=1; end end end %Calculate the values for the positions of plate 2 using the following %options with lsqnonlin options=optimset('TolFun',1E-18,'TolX',1E- 10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac obPattern',jac2); [XX,resnorm,residual,exitflag,output,lambda,jacobian2]=lsqnonlin(anon 2,X00,[],[],options); %Calculate the differences between the true and calculated values for midx, %midy, midz, alpha, beta and gamma difference=XX-truevalues; midx_difference=difference(1:6:end); midy_difference=difference(2:6:end); midz_difference=difference(3:6:end); alpha_difference=difference(4:6:end); beta_difference=difference(5:6:end); gamma_difference=difference(6:6:end); Page 48 Patrick Gloster %Calculate the standard deviation of these differences midx_SD=std(midx_difference); midy_SD=std(midy_difference); midz_SD=std(midz_difference); alpha_SD=std(alpha_difference); beta_SD=std(beta_difference); gamma_SD=std(gamma_difference); %Calculate the mean of these differences midx_mean=mean(midx_difference); midy_mean=mean(midy_difference); midz_mean=mean(midz_difference); alpha_mean=mean(alpha_difference); beta_mean=mean(beta_difference); gamma_mean=mean(gamma_difference); %Plot histograms for the results.. Histograms show the differences between the true and %calculated values of midx, midy, midz, alpha, beta and gamma h1=figure('Name','calibration bowtie 3 results','NumberTitle','off','Units'... ,'normalized','Position',[0 0 1 0.8]); subplot(2,3,1); hist(midx_difference,20) title('midx difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',midx_mean,midx_SD),'Location','SouthOutside') subplot(2,3,2); hist(midy_difference,20) title('midy difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',midy_mean,midy_SD),'Location','SouthOutside') subplot(2,3,3); hist(midz_difference,20) title('midz difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',midz_mean,midz_SD),'Location','SouthOutside') subplot(2,3,4); hist(alpha_difference,20) title('alpha difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',alpha_mean,alpha_SD),'Location','SouthOutside') subplot(2,3,5); hist(beta_difference,20) title('beta difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',beta_mean,beta_SD),'Location','SouthOutside') subplot(2,3,6); hist(gamma_difference,20) title('gamma difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',gamma_mean,gamma_SD),'Location','SouthOutside') Page 49 Patrick Gloster CALIBRATIONBOWTIE4.M function calibrationbowtie4(N,M,repeats) format long midx_SD=zeros(repeats,1); midy_SD=zeros(repeats,1); midz_SD=zeros(repeats,1); alpha_SD=zeros(repeats,1); beta_SD=zeros(repeats,1); gamma_SD=zeros(repeats,1); midx_mean=zeros(repeats,1); midy_mean=zeros(repeats,1); midz_mean=zeros(repeats,1); alpha_mean=zeros(repeats,1); beta_mean=zeros(repeats,1); gamma_mean=zeros(repeats,1); %This sets the seed for the random number by the clock time=clock; seed=time(6)+time(5)+time(4)+time(3)+time(2)+time(1); rand('seed',seed); randn('seed',seed); %These are the distances between retroreflectors on plate 2 d1=0.1; d2=0.1; initialerror=1E-3; %This is the error on our initial estimates of the positions, 1mm lasererror=1E-9; %This is the error on the laser measurements, 1nm distancebetweenplates=1; %This is the rough distance between the plates along the Z axis, 1 metre angleerror=(pi/180); %This is the error on how well we can determine the angle, currently 1 degree %These are the ranges of movement along each axis that plate 2 has platerangeinX=0.01; platerangeinY=0.01; Page 50 Patrick Gloster platerangeinZ=0.01; %These are the sizes of the angles that plate 2 can rotate through alphavariation=2.*pi/9; betavariation=2.*pi/9; gammavariation=2.*pi/9; %Define the points on plate 1 a1=0; b1=0; c1=0; a2=0; b2=0.02; c2=0; a3=0; b3=0.08; c3=0; a4=0; b4=0.1; c4=0; a5=-0.05; b5=0.045; c5=0; a6=-0.03; b6=0.045; c6=0; a7=0.03; b7=0.05; c7=0; a8=0.05; b8=0.05; c8=0; for count=1:repeats %Generate the midpoint values of plate 2 for N positions midx=platerangeinX.*2.*(rand(N,1)-0.5); midy=(b4/2)+platerangeinY.*2.*(rand(N,1)-0.5); midz=distancebetweenplates+platerangeinZ.*2.*(rand(N,1)-0.5); %Generate the angles of orientation of plate 2 for each of N positions alpha=2.*(rand(N,1)-0.5).*alphavariation; beta=2.*(rand(N,1)-0.5).*betavariation; gamma=2.*(rand(N,1)-0.5).*gammavariation; %Generate the true positions of points 1, 2, 3 and 4 x1=midx+(d1/2).*cos(alpha).*sin(gamma); Page 51 Patrick Gloster y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); truevalues=[midx,midy,midz,alpha,beta,gamma]; truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4]; truevalues=truevalues'; truevalues=truevalues(:); truepositions=truepositions'; truepositions=truepositions(:); %Generate our starting approximations for a's, b's and c's. Remember, a1=0, %b1=0, c1=0, a4=0, c4=0 and c6=0 so aren't included as variables approx_a2=a2+randn.*initialerror; approx_b2=b2+randn.*initialerror; approx_c2=c2+randn.*initialerror; approx_a3=a3+randn.*initialerror; approx_b3=b3+randn.*initialerror; approx_c3=c3+randn.*initialerror; approx_b4=b4+randn.*initialerror; approx_a5=a5+randn.*initialerror; approx_b5=b5+randn.*initialerror; approx_a6=a6+randn.*initialerror; approx_b6=b6+randn.*initialerror; approx_c6=c6+randn.*initialerror; approx_a7=a7+randn.*initialerror; approx_b7=b7+randn.*initialerror; approx_c7=c7+randn.*initialerror; approx_a8=a8+randn.*initialerror; approx_b8=b8+randn.*initialerror; approx_c8=c8+randn.*initialerror; Page 52 Patrick Gloster %Generate our initial approximations to d1 and d2: approx_d1=d1+randn.*initialerror; approx_d2=d2+randn.*initialerror; %Generate our approximations to the midpoints of plate 2 approx_midx=midx+randn(N,1).*initialerror; approx_midy=midy+randn(N,1).*initialerror; approx_midz=midz+randn(N,1).*initialerror; %Generate our approximations to the angles of orientation of plate 2 approx_alpha=alpha+(randn(N,1).*angleerror); approx_beta=beta+(rand(N,1).*angleerror); approx_gamma=gamma+(rand(N,1).*angleerror); approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro x_beta,approx_gamma]; approx_values=approx_values'; approx_values=approx_values(:); %Calculate the true distances D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); D4=sqrt((x2).^2+(y2).^2+(z2).^2); D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); %Generate the measurements by adding a random error to each distance L1=D1+randn(N,1).*lasererror; L2=D2+randn(N,1).*lasererror; L3=D3+randn(N,1).*lasererror; L4=D4+randn(N,1).*lasererror; L5=D5+randn(N,1).*lasererror; L6=D6+randn(N,1).*lasererror; L7=D7+randn(N,1).*lasererror; L8=D8+randn(N,1).*lasererror; measurements=[L1,L2,L3,L4,L5,L6,L7,L8]; %Form a vector of our initial estimates of all variables X0=[approx_a2;approx_b2;approx_c2;approx_a3;approx_b3;approx_c3;appro x_b4;... approx_a5;approx_b5;approx_a6;approx_b6;approx_c6;approx_a7;approx_b7 ;... approx_c7;approx_a8;approx_b8;approx_c8;approx_d1;approx_d2;approx_va lues]; %We need to use an anonymous function as the objective function because %superhappyfun has 2 input arguments and lsqnonlin can only minimize Page 53 Patrick Gloster %functions with 1 input argument anon = @(X) superhappyfun(measurements,X); %Set lsqnonlin to do one iteration. This will quickly perform the %minimzation with just one iteration, giving us a result with very bad %accuracy. However, what we are interested in is the jacobian which is %returned. Using this jacobian the sparsity structure is obtained and then %handed to lsqnonlin for the next minimization, which has options suitable %for our minimization and will run much faster using the sparsity %structure. options=optimset('MaxIter',1); [dummy_X,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,... dummy_lambda,dummy_jacobian]=lsqnonlin(anon,X0,[],[],options); %Derive the sparsity structure jac=dummy_jacobian; for a=1:8.*N for b=1:6.*N+20 if jac(a,b)~=0 jac(a,b)=1; end end end %Set the options for lsqnonlin suitable for a good minimization options=optimset('TolFun',1E-18,'TolX',1E- 10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac obPattern',jac); %Use the minimization function lsqnonlin to calculate values for %the positions of the launch heads on plate 1. These values are only %calculated once [X,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(anon,X 0,[],[],options); calibrations=X(1:20); clear midx midy midz alpha beta gamma x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4 clear truevalues truepositions approx_midx approx_midy approx_midz clear approx_alpha approx_beta approx_gamma approx_values clear D1 D2 D3 D4 D5 D6 D7 D8 L1 L2 L3 L4 L5 L6 L7 L8 measurements %Now generate M plate 2 positions and, using the calibrations we have just %worked out, calculate their positions %Generate the midpoint values of plate 2 for M positions Page 54 Patrick Gloster midx=platerangeinX.*2.*(rand(M,1)-0.5); midy=(b4/2)+platerangeinY.*2.*(rand(M,1)-0.5); midz=distancebetweenplates+platerangeinZ.*2.*(rand(M,1)-0.5); %Generate the angles of orientation of plate 2 for each of M positions alpha=2.*(rand(M,1)-0.5).*alphavariation; beta=2.*(rand(M,1)-0.5).*betavariation; gamma=2.*(rand(M,1)-0.5).*gammavariation; %Generate the true positions of points 1, 2, 3 and 4 x1=midx+(d1/2).*cos(alpha).*sin(gamma); y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); truevalues=[midx,midy,midz,alpha,beta,gamma]; truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4]; truevalues=truevalues'; truevalues=truevalues(:); truepositions=truepositions'; truepositions=truepositions(:); %Generate our approximations to the midpoints of plate 2 approx_midx=midx+randn(M,1).*initialerror; approx_midy=midy+randn(M,1).*initialerror; approx_midz=midz+randn(M,1).*initialerror; %Generate our approximations to the angles of orientation of plate 2 approx_alpha=alpha+(randn(M,1).*angleerror); approx_beta=beta+(rand(M,1).*angleerror); approx_gamma=gamma+(rand(M,1).*angleerror); approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro x_beta,approx_gamma]; approx_values=approx_values'; approx_values=approx_values(:); %Calculate the true distances D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); Page 55 Patrick Gloster D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); D4=sqrt((x2).^2+(y2).^2+(z2).^2); D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); %Generate the measurements by adding a random error to each distance L1=D1+randn(M,1).*lasererror; L2=D2+randn(M,1).*lasererror; L3=D3+randn(M,1).*lasererror; L4=D4+randn(M,1).*lasererror; L5=D5+randn(M,1).*lasererror; L6=D6+randn(M,1).*lasererror; L7=D7+randn(M,1).*lasererror; L8=D8+randn(M,1).*lasererror; measurements=[L1,L2,L3,L4,L5,L6,L7,L8]; X00=approx_values; anon2 = @(XX) calibrationfun(measurements,XX,calibrations); %Obtain the sparsity structure for this next minimization options=optimset('MaxIter',1); [dummy_XX,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,.. . dummy_lambda,dummy_jacobian2]=lsqnonlin(anon2,X00,[],[],options); jac2=dummy_jacobian2; for a=1:8.*M for b=1:6.*M if jac2(a,b)~=0 jac2(a,b)=1; end end end %Calculate the values for the positions of plate 2 using the following %options with lsqnonlin options=optimset('TolFun',1E-18,'TolX',1E- 10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac obPattern',jac2); [XX,resnorm,residual,exitflag,output,lambda,jacobian2]=lsqnonlin(anon 2,X00,[],[],options); %Calculate the differences between the true and calculated values for midx, %midy, midz, alpha, beta and gamma difference=XX-truevalues; midx_difference=difference(1:6:end); Page 56 Patrick Gloster midy_difference=difference(2:6:end); midz_difference=difference(3:6:end); alpha_difference=difference(4:6:end); beta_difference=difference(5:6:end); gamma_difference=difference(6:6:end); %Calculate the standard deviation of these differences midx_SD(count)=std(midx_difference); midy_SD(count)=std(midy_difference); midz_SD(count)=std(midz_difference); alpha_SD(count)=std(alpha_difference); beta_SD(count)=std(beta_difference); gamma_SD(count)=std(gamma_difference); %Calculate the mean of these differences midx_mean(count)=mean(midx_difference); midy_mean(count)=mean(midy_difference); midz_mean(count)=mean(midz_difference); alpha_mean(count)=mean(alpha_difference); beta_mean(count)=mean(beta_difference); gamma_mean(count)=mean(gamma_difference); calibration_numer=count end mean_midx_mean=mean(midx_mean); mean_midy_mean=mean(midy_mean); mean_midz_mean=mean(midz_mean); mean_alpha_mean=mean(alpha_mean); mean_beta_mean=mean(beta_mean); mean_gamma_mean=mean(gamma_mean); SD_midx_mean=std(midx_mean); SD_midy_mean=std(midy_mean); SD_midz_mean=std(midz_mean); SD_alpha_mean=std(alpha_mean); SD_beta_mean=std(beta_mean); SD_gamma_mean=std(gamma_mean); SD_midx_SD=std(midx_SD); SD_midy_SD=std(midy_SD); SD_midz_SD=std(midz_SD); SD_alpha_SD=std(alpha_SD); SD_beta_SD=std(beta_SD); SD_gamma_SD=std(gamma_SD); mean_midx_SD=mean(midx_SD); mean_midy_SD=mean(midy_SD); mean_midz_SD=mean(midz_SD); mean_alpha_SD=mean(alpha_SD); mean_beta_SD=mean(beta_SD); mean_gamma_SD=mean(gamma_SD); %Plot histograms for the results.. Histograms show the differences between the true and %calculated values of midx, midy, midz, alpha, beta and gamma Page 57 Patrick Gloster h1=figure('Name','means','NumberTitle','off','Units'... ,'normalized','Position',[0 0 1 0.8]); subplot(2,3,1); hist(midx_mean,20) title('midx mean'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_midx_mean,SD_midx_mean),'Location','SouthOutside') subplot(2,3,2); hist(midy_mean,20) title('midy mean'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_midy_mean,SD_midy_mean),'Location','SouthOutside') subplot(2,3,3); hist(midz_mean,20) title('midz mean'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_midz_mean,SD_midz_mean),'Location','SouthOutside') subplot(2,3,4); hist(alpha_mean,20) title('alpha mean'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_alpha_mean,SD_alpha_mean),'Location','SouthOutside ') subplot(2,3,5); hist(beta_mean,20) title('beta mean'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_beta_mean,SD_beta_mean),'Location','SouthOutside') subplot(2,3,6); hist(gamma_mean,20) title('gamma mean'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_gamma_mean,SD_gamma_mean),'Location','SouthOutside ') h2=figure('Name','Standard deviations','NumberTitle','off','Units'... ,'normalized','Position',[0 0 1 0.8]); subplot(2,3,1); hist(midx_SD,20) title('midx SD'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_midx_SD,SD_midx_SD),'Location','SouthOutside') subplot(2,3,2); hist(midy_SD,20) title('midy SD'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_midy_SD,SD_midy_SD),'Location','SouthOutside') subplot(2,3,3); hist(midz_SD,20) title('midz SD'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_midz_SD,SD_midz_SD),'Location','SouthOutside') Page 58 Patrick Gloster subplot(2,3,4); hist(alpha_SD,20) title('alpha SD'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_alpha_SD,SD_alpha_SD),'Location','SouthOutside') subplot(2,3,5); hist(beta_SD,20) title('beta SD'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_beta_SD,SD_beta_SD),'Location','SouthOutside') subplot(2,3,6); hist(gamma_SD,20) title('gamma SD'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_gamma_SD,SD_gamma_SD),'Location','SouthOutside') Page 59 Patrick Gloster SUPERHAPPYFUN.M – The objective function for calibrationbowtie3 and 4, and calibrations function out=superhappyfun(measurements,X) %This is the objective function for threeDbowtie. %First we define the current values of the a's, b's and c's on plate 1 from %the vector X: a2=X(1); b2=X(2); c2=X(3); a3=X(4); b3=X(5); c3=X(6); b4=X(7); a5=X(8); b5=X(9); a6=X(10); b6=X(11); c6=X(12); a7=X(13); b7=X(14); c7=X(15); a8=X(16); b8=X(17); c8=X(18); %Define the current values of d1 and d2 from the vector X: d1=X(19); d2=X(20); %Define the current values of midx, midy, midz, alpha, beta and gamma from %the vector X: midx=X(21:6:end); midy=X(22:6:end); midz=X(23:6:end); alpha=X(24:6:end); beta=X(25:6:end); gamma=X(26:6:end); %Calculate the positions of points 1, 2, 3 and 4 based on our current %values for midx, midy, midz, alpha, beta and gamma: x1=midx+(d1/2).*cos(alpha).*sin(gamma); y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); Page 60 Patrick Gloster x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); %Calculate what measurements we would observe if the 4 points were at these %calculated positions. These are the 'predicted measurements': M1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); M2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); M3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); M4=sqrt((x2).^2+(y2).^2+(z2).^2); M5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); M6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); M7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); M8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); predicted=[M1,M2,M3,M4,M5,M6,M7,M8]; %Calculate the difference between these 'predicted' measurements and the %true measurements: out=predicted-measurements; Page 61 Patrick Gloster CALIBRATIONFUN.M – A second objective function used by calbration3 and 4 after they have calibrated plate 1 function out=calibrationfun(measurements,X,calibrations) %This is the objective function for threeDbowtie. %First we define the current values of the a's, b's and c's on plate 1 from %the vector calibrations: a2=calibrations(1); b2=calibrations(2); c2=calibrations(3); a3=calibrations(4); b3=calibrations(5); c3=calibrations(6); b4=calibrations(7); a5=calibrations(8); b5=calibrations(9); a6=calibrations(10); b6=calibrations(11); c6=calibrations(12); a7=calibrations(13); b7=calibrations(14); c7=calibrations(15); a8=calibrations(16); b8=calibrations(17); c8=calibrations(18); %Define the current values of d1 and d2 from the vector calibrations: d1=calibrations(19); d2=calibrations(20); %Define the current values of midx, midy, midz, alpha, beta and gamma from %the vector X: midx=X(1:6:end); midy=X(2:6:end); midz=X(3:6:end); alpha=X(4:6:end); beta=X(5:6:end); gamma=X(6:6:end); %Calculate the positions of points 1, 2, 3 and 4 based on our current %values for midx, midy, midz, alpha, beta and gamma: x1=midx+(d1/2).*cos(alpha).*sin(gamma); y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); Page 62 Patrick Gloster x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); %Calculate what measurements we would observe if the 4 points were at these %calculated positions. These are the 'predicted measurements': M1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); M2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); M3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); M4=sqrt((x2).^2+(y2).^2+(z2).^2); M5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); M6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); M7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); M8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); predicted=[M1,M2,M3,M4,M5,M6,M7,M8]; %Calculate the difference between these 'predicted' measurements and the %true measurements: out=predicted-measurements; Page 63 Patrick Gloster CALIBRATIONS.M function calibrations(N) format long %This sets the seed for the random number by the clock time=clock; seed=time(6)+time(5)+time(4)+time(3)+time(2)+time(1); rand('seed',seed); randn('seed',seed); %These are the distances between retroreflectors on plate 2 d1=0.1; d2=0.1; initialerror=1E-3; %This is the error on our initial estimates of the positions, 1mm lasererror=1E-9; %This is the error on the laser measurements, 1nm distancebetweenplates=1; %This is the rough distance between the plates along the Z axis, 1 metre angleerror=(pi/180); %This is the error on how well we can determine the angle, currently 1 degree %These are the ranges of movement along each axis that plate 2 has platerangeinX=0.01; platerangeinY=0.01; platerangeinZ=0.01; %These are the sizes of the angles that plate 2 can rotate through alphavariation=2.*pi/9; betavariation=2.*pi/9; gammavariation=2.*pi/9; %Define the points on plate 1 a1=0; b1=0; c1=0; a2=0; b2=0.02; c2=0; a3=0; b3=0.08; c3=0; a4=0; b4=0.1; c4=0; a5=-0.05; b5=0.045; c5=0; a6=-0.03; b6=0.045; Page 64 Patrick Gloster c6=0; a7=0.03; b7=0.05; c7=0; a8=0.05; b8=0.05; c8=0; for count=1:50 %Generate the midpoint values of plate 2 for N positions midx=platerangeinX.*2.*(rand(N,1)-0.5); midy=(b4/2)+platerangeinY.*2.*(rand(N,1)-0.5); midz=distancebetweenplates+platerangeinZ.*2.*(rand(N,1)-0.5); %Generate the angles of orientation of plate 2 for each of N positions alpha=2.*(rand(N,1)-0.5).*alphavariation; beta=2.*(rand(N,1)-0.5).*betavariation; gamma=2.*(rand(N,1)-0.5).*gammavariation; %Generate the true positions of points 1, 2, 3 and 4 x1=midx+(d1/2).*cos(alpha).*sin(gamma); y1=midy-(d1/2).*cos(alpha).*cos(gamma); z1=midz+(d1/2).*sin(alpha); x2=midx-(d1/2).*cos(alpha).*sin(gamma); y2=midy+(d1/2).*cos(alpha).*cos(gamma); z2=midz-(d1/2).*sin(alpha); x3=midx+(d2/2).*cos(beta).*cos(gamma); y3=midy+(d2/2).*cos(beta).*sin(gamma); z3=midz-(d2/2).*sin(beta); x4=midx-(d2/2).*cos(beta).*cos(gamma); y4=midy-(d2/2).*cos(beta).*sin(gamma); z4=midz+(d2/2).*sin(beta); truevalues=[midx,midy,midz,alpha,beta,gamma]; truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4]; truevalues=truevalues'; truevalues=truevalues(:); truepositions=truepositions'; truepositions=truepositions(:); %Generate our starting approximations for a's, b's and c's. Remember, a1=0, %b1=0, c1=0, a4=0, c4=0 and c6=0 so aren't included as variables Page 65 Patrick Gloster approx_a2=a2+randn.*initialerror; approx_b2=b2+randn.*initialerror; approx_c2=c2+randn.*initialerror; approx_a3=a3+randn.*initialerror; approx_b3=b3+randn.*initialerror; approx_c3=c3+randn.*initialerror; approx_b4=b4+randn.*initialerror; approx_a5=a5+randn.*initialerror; approx_b5=b5+randn.*initialerror; approx_a6=a6+randn.*initialerror; approx_b6=b6+randn.*initialerror; approx_c6=c6+randn.*initialerror; approx_a7=a7+randn.*initialerror; approx_b7=b7+randn.*initialerror; approx_c7=c7+randn.*initialerror; approx_a8=a8+randn.*initialerror; approx_b8=b8+randn.*initialerror; approx_c8=c8+randn.*initialerror; %Generate our initial approximations to d1 and d2: approx_d1=d1+randn.*initialerror; approx_d2=d2+randn.*initialerror; %Generate our approximations to the midpoints of plate 2 approx_midx=midx+randn(N,1).*initialerror; approx_midy=midy+randn(N,1).*initialerror; approx_midz=midz+randn(N,1).*initialerror; %Generate our approximations to the angles of orientation of plate 2 approx_alpha=alpha+(randn(N,1).*angleerror); approx_beta=beta+(rand(N,1).*angleerror); approx_gamma=gamma+(rand(N,1).*angleerror); approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro x_beta,approx_gamma]; approx_values=approx_values'; approx_values=approx_values(:); %Calculate the true distances D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2); D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2); D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2); D4=sqrt((x2).^2+(y2).^2+(z2).^2); D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2); D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2); D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2); D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2); Page 66 Patrick Gloster %Generate the measurements by adding a random error to each distance L1=D1+randn(N,1).*lasererror; L2=D2+randn(N,1).*lasererror; L3=D3+randn(N,1).*lasererror; L4=D4+randn(N,1).*lasererror; L5=D5+randn(N,1).*lasererror; L6=D6+randn(N,1).*lasererror; L7=D7+randn(N,1).*lasererror; L8=D8+randn(N,1).*lasererror; measurements=[L1,L2,L3,L4,L5,L6,L7,L8]; %Form a vector of our initial estimates of all variables X0=[approx_a2;approx_b2;approx_c2;approx_a3;approx_b3;approx_c3;appro x_b4;... approx_a5;approx_b5;approx_a6;approx_b6;approx_c6;approx_a7;approx_b7 ;... approx_c7;approx_a8;approx_b8;approx_c8;approx_d1;approx_d2;approx_va lues]; %We need to use an anonymous function as the objective function because %superhappyfun has 2 input arguments and lsqnonlin can only minimize %functions with 1 input argument anon = @(X) superhappyfun(measurements,X); %Set lsqnonlin to do one iteration. This will quickly perform the %minimzation with just one iteration, giving us a result with very bad %accuracy. However, what we are interested in is the jacobian which is %returned. Using this jacobian the sparsity structure is obtained and then %handed to lsqnonlin for the next minimization, which has options suitable %for our minimization and will run much faster using the sparsity %structure. options=optimset('MaxIter',1); [dummy_X,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,... dummy_lambda,dummy_jacobian]=lsqnonlin(anon,X0,[],[],options); %Derive the sparsity structure jac=dummy_jacobian; for a=1:8.*N for b=1:6.*N+20 if jac(a,b)~=0 jac(a,b)=1; end end end Page 67 Patrick Gloster %Set the options for lsqnonlin suitable for a good minimization options=optimset('TolFun',1E-18,'TolX',1E- 10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac obPattern',jac); %Use the minimization function lsqnonlin to calculate values for %the positions of the launch heads on plate 1. These values are only %calculated once [X,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(anon,X 0,[],[],options); calibrations=X(1:20); a2_difference(count)=a2-calibrations(1); b2_difference(count)=b2-calibrations(2); c2_difference(count)=c2-calibrations(3); a3_difference(count)=a3-calibrations(4); b3_difference(count)=b3-calibrations(5); c3_difference(count)=c3-calibrations(6); b4_difference(count)=b4-calibrations(7); a5_difference(count)=a5-calibrations(8); b5_difference(count)=b5-calibrations(9); a6_difference(count)=a6-calibrations(10); b6_difference(count)=b6-calibrations(11); c6_difference(count)=c6-calibrations(12); a7_difference(count)=a7-calibrations(13); b7_difference(count)=b7-calibrations(14); c7_difference(count)=c7-calibrations(15); a8_difference(count)=a8-calibrations(16); b8_difference(count)=b8-calibrations(17); c8_difference(count)=c8-calibrations(18); d1_difference(count)=d1-calibrations(19); d2_difference(count)=d2-calibrations(20); repeat=count end mean_a2=mean(a2_difference); SD_a2=std(a2_difference); mean_b2=mean(b2_difference); SD_b2=std(b2_difference); mean_c2=mean(c2_difference); SD_c2=std(c2_difference); mean_a3=mean(a3_difference); SD_a3=std(a3_difference); mean_b3=mean(b3_difference); SD_b3=std(b3_difference); mean_c3=mean(c3_difference); SD_c3=std(c3_difference); mean_b4=mean(b4_difference); SD_b4=std(b4_difference); mean_a5=mean(a5_difference); SD_a5=std(a5_difference); mean_b5=mean(b5_difference); SD_b5=std(b5_difference); mean_a6=mean(a6_difference); SD_a6=std(a6_difference); Page 68 Patrick Gloster mean_b6=mean(b6_difference); SD_b6=std(b6_difference); mean_c6=mean(c6_difference); SD_c6=std(c6_difference); mean_a7=mean(a7_difference); SD_a7=std(a7_difference); mean_b7=mean(b7_difference); SD_b7=std(b7_difference); mean_c7=mean(c7_difference); SD_c7=std(c7_difference); mean_a8=mean(a8_difference); SD_a8=std(a8_difference); mean_b8=mean(b8_difference); SD_b8=std(b8_difference); mean_c8=mean(c8_difference); SD_c8=std(c8_difference); mean_d1=mean(d1_difference); SD_d1=std(d1_difference); mean_d2=mean(d2_difference); SD_d2=std(d2_difference); h2=figure('Name','calibrations - a differences','NumberTitle','off','Units'... ,'normalized','Position',[0 0 1 0.8]); subplot(2,3,1); hist(a2_difference,5) title('a2 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_a2,SD_a2),'Location','SouthOutside') subplot(2,3,2); hist(a3_difference,5) title('a3 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_a3,SD_a3),'Location','SouthOutside') subplot(2,3,3); hist(a5_difference,5) title('a5 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_a5,SD_a5),'Location','SouthOutside') subplot(2,3,4); hist(a6_difference,5) title('a6 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_a6,SD_a6),'Location','SouthOutside') subplot(2,3,5); hist(a7_difference,5) title('a7 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_a7,SD_a7),'Location','SouthOutside') subplot(2,3,6); hist(a8_difference,5) title('a8 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_a8,SD_a8),'Location','SouthOutside') Page 69 Patrick Gloster %Plot histograms of the differences between the true and calculated values %of b2, b3, b4, b5, b6, b7 and b8 h3=figure('Name','calibrations - b differences','NumberTitle','off','Units'... ,'normalized','Position',[0 0 1 0.8]); subplot(2,4,1); hist(b2_difference,5) title('b2 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b2,SD_b2),'Location','SouthOutside') subplot(2,4,2); hist(b3_difference,5) title('b3 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b3,SD_b3),'Location','SouthOutside') subplot(2,4,3); hist(b4_difference,5) title('b4 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b4,SD_b4),'Location','SouthOutside') subplot(2,4,4); hist(b5_difference,5) title('b5 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b5,SD_b5),'Location','SouthOutside') subplot(2,4,5); hist(b6_difference,5) title('b6 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b6,SD_b6),'Location','SouthOutside') subplot(2,4,6); hist(b7_difference,5) title('b7 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b7,SD_b7),'Location','SouthOutside') subplot(2,4,7); hist(b8_difference,5) title('b8 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_b8,SD_b8),'Location','SouthOutside') %Plot histograms of the differences between the true and calculated values %of c2, c3, c6, c7 and c8 h4=figure('Name','calibrations - c differences','NumberTitle','off','Units'... ,'normalized','Position',[0 0 1 0.8]); subplot(2,3,1); Page 70 Patrick Gloster hist(c2_difference,5) title('c2 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_c2,SD_c2),'Location','SouthOutside') subplot(2,3,2); hist(c3_difference,5) title('c3 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_c3,SD_c3),'Location','SouthOutside') subplot(2,3,3); hist(c6_difference,5) title('c6 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_c6,SD_c6),'Location','SouthOutside') subplot(2,3,4); hist(c7_difference,5) title('c7 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_c7,SD_c7),'Location','SouthOutside') subplot(2,3,5); hist(c8_difference,5) title('c8 difference'); legend(sprintf('Mean=%e\nStandard deviation=%e',mean_c8,SD_c8),'Location','SouthOutside') Page 71