If the field of the entire machine obeys B(x,y,z+a) = e^(ka)B(x,y,z), then the tunes ought to be fixed, but this generally won't be true far out in the fringe fields in any practical coil shape.

Exactly, so moving from constant m value to constant tunes optimisation is the only practical way here

In fact I think the only magnetic fields that scale perfectly everywhere also increase exponentially forever upwards, so are only of theoretical interest

Is your magnet just bare coils? If so you might want to consider an iron collar like in a traditional cos-theta superconducting magnet. "Collar" is not really the right word here with this structure, but the idea is to have iron surrounding the the entire coil structure. I believe it's made thick enough that it returns the entirety of the flux. This will help reduce the influence of the surrounding magnets. The problem is that the iron will probably screw up your fields-from-Biot-Savart method (I think it basically works for SC magnets because the iron near the coils is saturated). You don't avoid the magnet-to-magnet influence completely, but with that large of a range of influence, I think it would be worth it. It will also shield you from being influenced by the surroundings.

yes, only coils are considered at the moment, as analytical solution can be found in this way. I agree that shielding will change the field profile and will thus have to be reshaped with Opera. We are still in the feasibility study, and the winding method is critical at this point to confirm the possibility to manufacture these magnets

This is actually in our plan, there'll eventually be an iterative design process (hopefully automated) where the winding density is retuned to take account of the iron shielding.  This is why we're starting to use computing cluster time for OPERA simulations...