Not looking to spend a lot of run time awaiting an answer, I ran

clc,clearvars

global phi availOps stack availDigs bestForm bestVal form

phi=(1+sqrt(5))/2;

digSet=1:9;

availOps='+-*/';

stack=[];

availDigs=digSet;

bestVal=9999999;

form='';

findIt()

function findIt()

global phi availOps stack availDigs bestForm bestVal form 

for whDig=1:length(availDigs)

   dig=availDigs(whDig);

   availDigs(whDig)=[];

   stack=[stack dig];

   form=[form char(string(dig))];

   findIt();

   stack(length(stack))=[];

   form(length(form))=[];

   availDigs=[availDigs(1:whDig-1) dig availDigs(whDig:end)];

end

if length(stack)>1

  for whOp=1:length(availOps)

    Op=availOps(whOp);

    saveform=form;

    form=[form Op];

    if Op~='/' || stack(end)~=0

      v=eval([num2str(stack(end-1),16) Op num2str(stack(end),16)]);

      if abs(phi-v)<abs(phi-bestVal)

         bestVal=v;

         bestForm=form(length(stack)-1:end); % before stack reduction

         disp([v phi-v])

         disp(['     ' bestForm])

         disp(' ')

      end

      stacksave=stack;

      stack=[stack(1:end-2) v];

      if ~isempty(availDigs)

        findIt();

      end

      stack=stacksave;     

    end

    form=saveform;

  end

end

end

for a few minutes. It produces the RPN (Reverse Polish Notation) formula for the computation. It produces some results using only part of the set of 9 digits (e.g., (7+8)/9 was the algebraic equivalent of an early one), but the later, more accurate, ones used all nine.  The last one produced before I terminated the program gave the following:

          The approximation         How far off
          1.61801242236025      2.15663896463436e-05
          
      The RPN formula     
     1234569*8-/+*--7/

     

Converted to Algebraic form, that's 

(1-(2-(3*(4+5/(6*9-8)))))/7.

I also tried reversing the order of the attempted digits, and using random order several times, but this remained the best I found.