If you for whatever reason cannot use one of the existing solver implementations, you can improve the estimation heuristic by not just counting how many cells are incorrectly positioned but rather summing the distance(whatever that might be) from their correct positions.
0. State representation
You can represent position of a cell by 3D vector(=array) of 3 values in range <-1,1> that is from -1 left/bottom, to 0 for center to right/top for +1. For example the [1,1,0] is for upper-right cell of front face. Center of cube is the origin. Note: face center cells are fixed so the position is orientation invariant as long as you always orient face centers the same way.
1. Manhattan distance
One of the simplest metrics is so called Manhattan distance, where you just add absolute differences of corresponding elements. For example manhattan distance from int[] v0 = [0,1,-1] to int[] v1 = int [1,1,1] is
var manhattan = abs(v0[0]-v1[0]) + abs(v0[1]-v1[1]) + abs(v0[2]-v1[2])
the solution level estimate is simply sum of Manhattan distance to correct position for each cell. While this is likely to be big improvement, it still does not take in account that the cube wraps around.
If I am not mistaken the Manhattan distance was used as heuristic for IDA*, one of the first optimal solvers (in terms of least turns).
2. Swap & negate
As you might noticed, the actual distance to solution is corresponding to number (quarter)turns rather than physical distance on the cube. Instead, we should perform turns themselves - 3D rotation of the cell position vectors around axes which would lead to matrix multiplication...
...However as it turns out our case can be greatly simplified. Thanks to choosing <-1,1> instead of <0,2> we can use the well known swap & negate "trick", a way calculate perpendicular vectors. Fixing one element of vector, swapping the other two and negating one(which one controls the direction) of those will perform the quarter turn on face perpendicular to the fixed element. The procedure could look like this.
//calculates position of cell after a quarter-turn around given axis
void Rotate3DAligned(int[] vec, int axis, bool ccw)
{
var first = vec[(i-1)%3];
var second = vec[(i+1)%3];
if(ccw) first = -first;
else second = -second;
vec[(i-1)%3] = second;
vec[(i+1)%3] = first ;
}
Taking the current position of each cell and the correct one you need to determine how many turns (rotations - swap&negates) you need to take in order get it into the correct postion. This can be trivially solved by bruteforce since the max number is rather low, or better algorithmicaly - how many signs you need swap and how many values to "move" around.
2.1. Hash it
Count only distinct turn series i.e. save the series of turns into a HashSet or similar. This will allow detect many cases of "batch" moving cells, solving the one-turn-from solution problem.
disclaimer: algorithm not tested by implementation, code not checked
