Surface fitting with multiple 2d gaussian functions

Consider the functions:

f1=A1.*exp(-((x-x01).^2./(2.*sigmax1^2)+(y-y01).^2./(2.*sigmay1^2)))

f2=A2.*exp(-((x-x02).^2./(2.*sigmax2^2)+(y-y02).^2./(2.*sigmay2^2)))

f3=A3.*exp(-((x-x03).^2./(2.*sigmax3^2)+(y-y03).^2./(2.*sigmay3^2)))

with respective volumes given by

V1=2*pi*A1*sigmax1*sigmay1

V2=2*pi*A2*sigmax2*sigmay2

V3=2*pi*A3*sigmax3*sigmay3

Problem Statement:

Known:

1. Sum of f1+f2+f3 is known across domain [x,y]=meshgrid();
2. Ok to assume that f1+f2+f3=0 at the domain boundary
3. The volumes for V1, V2, V3 are known

Objective:

Find parameters A1,x01,y01,sigmax1,sigmay1,A2,x02,y02,sigmax2,sigmay2,A3,x03,y03,sigmax3,sigmay3

Any suggestion on how to solve such a problem would be highly appreciated.