Блог

Я пытаюсь определить гидравлический параметр (высота водной поверхности, WSE), который приводит к тому, что мой расчетный расход потока соответствует измеренному расходу в данном месте. У меня есть много таких мест для повторения, поэтому идеально создать цикл для запуска всех моих данных. На данный момент я просто тестирую это в заданном месте поперечного сечения.

То, что я хочу, чтобы optim() сделал, — это нашел наилучший WSE (высота водной поверхности), который входит в другие гидравлические уравнения (смоченная площадь, смоченный периметр) для оценки расхода потока. Прямо сейчас кажется, что optim() просто дает мне значение разрядки. Предполагая, что у меня просто неправильный порядок операций или что я указываю optim() не на ту вещь? Любое понимание очень ценится.

 # make up some dummy data: x lt;- seq(0, 50, length.out = 20) z lt;- c(seq(100,60,length.out = 8), 58,56,56,58, seq(60,100, length.out = 8)) plot(x,z) # v-shaped stream channel cross section  #remove NAs (lots of these in my real data): x lt;- x[!is.na(x)] z lt;- z[!is.na(z)]  # Mannings N: N lt;- 0.05  # Target Q (discharge, in cubic feet per second):  Q lt;- 400  # Slope (ft/ft) slope lt;- 0.004  # Create a function to minimize the difference between estimated and calculated Q: Q.calc lt;- function(WSE){    # calculate wetted area:  wet.area lt;- vector(length = length(x))  for(i in 1:(length(x)-1)){  if(max(z[i],z[i 1]) gt; WSE){ # assign 0 to locations above WSE  wet.area[i] lt;- 0  } else {  wet.area[i] lt;- (x[i 1] -x[i])*(WSE - (z[i 1]   z[i])/2)  }  }    wet.area lt;- sum(wet.area) # sum across channel    # calculate wetted perimeter:  wet.perim lt;- vector(length = length(z))  for(i in 1:(length(z)-1)){  if(max(z[i],z[i 1]) gt; WSE) { # assign 0 to locations above WSE  wet.perim[i] lt;- 0  } else { # calculate perimeter for locations below WSE  wet.perim[i] lt;- sqrt(((x[i 1]-x[i])^2) ((abs(z[i 1]-z[i])^2)))  }  }  wet.perim lt;- sum(wet.perim) # sum across channel    # hydraulic radius:  hyd.radius lt;- wet.area/wet.perim    # This is the equation we're actually 'optimizing':  # Calculating when Q.calc = Q observed:   (1.49/N)*wet.area*hyd.radius^(2/3)*slope^0.5 - Q # since we're minimizing we   #subtract the observed Q so that optimal Q.calc = Q.obs ??    }  optim(par = mean(z),   fn = Q.calc,   lower = min(z),   upper = max(z),  method = 'Brent',  control = list(maxit = 100000))  $par [1] 56  $value [1] -400  $counts function gradient   NA NA   $convergence [1] 0  $message NULL