Hi All,

I have made a number of post lately relating to this problem so I thought I would make a decent post and highlight exactly what my problem is so that hopefully I can get some feedback. Basically what I am doing is using the optimisation toolbox in matlab to find optimize the geometric potential of a system of nodes (put simply trying to get the nodes to disperse to find a minimum energy state).

The geometric potential of a system of uniformly weighted nodes is expressed by:

$\displaystyle G = \sum_{i=1}^{n}\frac{1}{d_{ij}^q}$

Where:

$\displaystyle n$ is the number of nodes in the system

$\displaystyle q$ is the power index

$\displaystyle d_{ij}$ is the distance between node $\displaystyle i$ and node $\displaystyle j$

Two forms of geometric potential are considered:

Global Geometric Potential $\displaystyle \left( G_G \right)$ and Relative Geometric Potential $\displaystyle \left( G_R \right)$

In the case of $\displaystyle G_G$, each node is considered to be connected to every other node. $\displaystyle G_R$ only considers the contributions from adjacent nodes that are connected directly. $\displaystyle G_R$ is the form used by the objective function in this example.

The figure below shows a set of 9 nodes placed at random inside a unit square. The black lines represent the connections between nodes and the blue lines represent the voronoi diagram of the points.

The problem is that the optimization routine is changing the values of $\displaystyle x_i$ and $\displaystyle y_i$ by too large a value per iteration. The connectivity as determined by the Delaunay triangulation is only valid inside each nodes Voronoi cell (in blue) so the maximum amount any node should move after only 1 iteration is to the edge of its respective Voronoi cell.

My understanding is that the optimization routine recognizes that moving node $\displaystyle 1$ and node $\displaystyle 4$ together to the point $\displaystyle \left(1,0\right)$ (for example) is a valid since the two nodes are not connected. The problem then arises that the Delaunay triangulation can no longer be recomputed since two nodes would now occupy the same position so there are duplicate nodes in the system.

I have tried writing a constraint function to limit how close 2 nodes can get to each other, but the step size taken by the routine appears to be the biggest issue that needs addressing. my problem in a nutshell is:

Is there a way to reduce the step size of the problem limit $\displaystyle x_i$ and $\displaystyle y_i$ to move by a maximum absolute value of say $\displaystyle \delta = 0.05$ per iteration?

The matlab code used is given here. Running the code without any parameters will use the default values shown in the figure above. I am useing MATLAB 2008b with the optimisation toolbox version 4.1.

Any advice on how to solve this issue would be greatly appreciated.Code:`function [xy,fval,eflag] = Relative_Geometric_Potential(xy_start,q)`

% Take the nodes xy_start and optimize their position inside a

% unit sqaure in order to minimize the relative geometric

% potential of the system subject to the power index q.

% > xy_start is a 1x2n vector where there are n nodes in the system.

% > q is a possitive integer.

if nargin < 2;q = 2;end%set q=2 by default.(Not that important for this example)

%Start position of nodes [x,y]

if nargin < 1 % use default

xy_start = [0.196 0.992 0.802 0.729 0.498 0.809 0.357 0.073 0.591,...%xvals

0.910 0.194 0.432 0.749 0.946 0.764 0.559 0.184 0.498];%yvals

end

n = length(xy_start)/2;%number of nodes

LBounds = zeros(1,2*n);%lower Bounds

UBounds = ones(1,2*n);%upper bounds

options = optimset('Display','iter','Algorithm','active-set');

f = @(xy)Gr(xy,q);%use power index of 2 for Geometric Potential Function

[xy,fval,eflag] = fmincon(f,xy_start,...

[],[],[],[],LBounds,UBounds,[],options);

xy = reshape(xy,length(xy)/2,2);%reshape to nx2 matrix

end

function fval = Gr(xy,q)

%xy is row vector = [x y]

%q is power index

xy = reshape(xy,length(xy)/2,2);%reshape to nx2 matrix

x = xy(:,1);y = xy(:,2);%pull column vectors

TRI = delaunay(x,y,{'QJ' 'Pp'});%determine connectivity between local nodes

tmp = [TRI(:,[1,2]);TRI(:,[2,3])];%find nodes i,j of individual element

tmp = sort(tmp,2);%sort so node i < node j

mem_ij = unique(tmp,'rows');%take unique rows to get connectivity

k = xy(mem_ij(:,1),:)-xy(mem_ij(:,2),:);%vector delta

dij = sqrt(sum(k.^2,2));%normal of k

foo = sum(dij.^(-q));%sum for final ans

fval = 2*foo;%account for dij -> dji connections (ie symetric dij matrix)

end

Kind Regards,

Elbarto