I would like to obtain an approximation of a function $1/(a x + b y)$ around $x=1$ in terms of a Laurent polynomial. I do not need the complete series; in fact perhaps the 1st order terms in $x$ and $x^{-1}$ may be sufficient for my purposes, leading to something of the form $\alpha_{-1} x^{-1} + \alpha_0 + \alpha_1 x$. However a Taylor expansion is not what I need. What is the correct way of obtaining such an expansion?
The next step is to expand in $y$ around 1 as well, obtaining something of the form $\alpha_{-1} x^{-1} + \beta_{-1} y^{-1} + \alpha_0 + \alpha_1 x + \beta_1 y$. Is this reasonable?
Expanding around $x=1$ and $y=1$ seem to be the wrong approaches to take here. The standard approach is to factor the larger (in magnitude) of $ax$ and $by$ out from the denominator, and use the geometric series summation formula in reverse (equivalently, the binomial expansion formula). So: $$\begin{align} \frac{1}{ax + by} & = \frac{1}{ax} \sum_{n=0}^\infty \left(-\frac{by}{ax}\right)^n \ \forall\, |ax| > |by|,\ \mathrm{and} \\ & = \frac{1}{by} \sum_{n=0}^\infty \left(-\frac{ay}{by}\right)^n \ \forall\, |ax| < |by| \end{align}.$$ If you want to expand about a different point, just shift your variables and apply the same formulae. For example: $$\begin{align} x' &= x - 1,\ \mathrm{and} \\ y' & = y + \frac{a}{b}. \end{align}$$
To more literally apply repeated expansions, as requested in the original post, you'll apply a binomial expansion to each term in one of the above formulae after applying another shift/change of variables. The formula for more general binomial expansions is: $$\begin{align}\frac{1}{(ax + b)^n} &= \frac{1}{(ax)^n} \sum_{m=0}^\infty \frac{\Gamma(-n+1)}{m! \Gamma(-n-m+1)} \left(\frac{b}{ax}\right)^m\ \forall\, |ax| > |b|,\ \mathrm{and} \\ & = \frac{1}{b^n} \sum_{m=0}^\infty \frac{\Gamma(-n+1)}{m! \Gamma(-n-m+1)} \left(\frac{ax}{b}\right)^m\ \forall\, |ax| < |b|. \end{align}$$