Subversion Repositories public iLand

Rev

Rev 1221 | Blame | Compare with Previous | Last modification | View Log | RSS feed

/********************************************************************************************
**    iLand - an individual based forest landscape and disturbance model
**    http://iland.boku.ac.at
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
**
**    This program is free software: you can redistribute it and/or modify
**    it under the terms of the GNU General Public License as published by
**    the Free Software Foundation, either version 3 of the License, or
**    (at your option) any later version.
**
**    This program is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**    GNU General Public License for more details.
**
**    You should have received a copy of the GNU General Public License
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************************************/

#include "abe_global.h"
#include "global.h"

#include "fmunit.h"

#include "forestmanagementengine.h"
#include "fmstand.h"
#include "scheduler.h"
#include "agent.h"
#include "agenttype.h"
#include "fomescript.h"
#include "fmtreelist.h"

namespace ABE {

/** @class FMUnit
    @ingroup abe
    The FMUnit class encapsulates a forest management unit, comprised of forest stands. Units are the base level at which
    the scheduling works.

  */



void FMUnit::aggregate()
{
    // loop over all stands
    // collect some data....
    double age=0.;
    double volume = 0.;
    double harvest = 0.;
    double totalarea = 0.;
    const QMultiMap<FMUnit*, FMStand*> &stands = ForestManagementEngine::instance()->stands();
    QMultiMap<FMUnit*, FMStand*>::const_iterator it = stands.constFind(this);
    while (it != stands.constEnd() && it.key()==this) {
        const FMStand *s = it.value();
        age += s->age() * s->area();
        volume += s->volume() * s->area();

        totalarea += s->area();
        ++it;
    }
    if (totalarea>0.) {
        age /= totalarea;
        volume /= totalarea;
        harvest /= totalarea;
    }
    qCDebug(abe) << "unit" << id() << "volume (m3/ha)" << volume << "age" << age << "planned harvest: todo";

}

QStringList FMUnit::info() const
{
    return QStringList() << QString("(accumulated) harvest: %1").arg(mRealizedHarvest)
                         << QString("MAI: %1").arg(mMAI)
                         << QString("HDZ: %1").arg(mHDZ)
                         << QString("average age: %1").arg(mMeanAge)
                         << QString("decadal plan: %1").arg(mAnnualHarvestTarget)
                         << QString("current plan: %1").arg(constScheduler()!=0?constScheduler()->harvestTarget():0.);

}

double FMUnit::annualThinningHarvest() const
{
    const QMultiMap<FMUnit*, FMStand*> &stands = ForestManagementEngine::instance()->stands();
    QMultiMap<FMUnit*, FMStand*>::const_iterator it = stands.constFind(const_cast<FMUnit*>(this));
    double harvested=0.;
    while (it != stands.constEnd() && it.key()==this) {
        FMStand *stand = it.value();
        harvested += stand->totalThinningHarvest();
        ++it;
    }
    return harvested;
}

FMUnit::FMUnit(const Agent *agent)
{
    mAgent = agent;
    mScheduler = 0;
    mAnnualHarvestTarget = -1.;
    mRealizedHarvest = 0.;
    mMAI = 0.; mHDZ = 0.; mMeanAge = 0.;
    mTotalArea = 0.; mTotalPlanDeviation = 0.;
    mTotalVolume = 0.;
    mAnnualHarvest = 0.;
    mNumberOfStands = 0;
    mU = 100, mThinningIntensityClass = 2, mSpeciesCompositionIndex = 0;
    mAverageMAI = 0.;


    //if (agent->type()->schedulerOptions().useScheduler)
    // explicit scheduler only for stands/units that include more than one stand
    mScheduler = new Scheduler(this);
}

FMUnit::~FMUnit()
{
    if (mScheduler)
        delete mScheduler;
}

void FMUnit::setId(const QString &id)
{
    mId = id;
}

void FMUnit::resetHarvestCounter()
{
    if (scheduler())
        scheduler()->resetHarvestCounter();
}

void FMUnit::managementPlanUpdate()
{
    const double period_length = 10.;
    // calculate the planned harvest in the next planning period (i.e., 10yrs).
    // this is the sum of planned operations that are already in the scheduler.
    double plan_thinning, plan_final;
    mScheduler->plannedHarvests(plan_final, plan_thinning);
    // the actual harvests of the last planning period
    //double realized = mRealizedHarvest;

    mRealizedHarvest = 0.; // reset
    mRealizedHarvestLastYear = 0.;


    // preparations:
    // MAI-calculation for all stands:
    double total_area = 0.;
    double age = 0.;
    double mai = 0.;
    double hdz = 0.;
    double volume = 0.;
    const QMultiMap<FMUnit*, FMStand*> &stands = ForestManagementEngine::instance()->stands();
    QMultiMap<FMUnit*, FMStand*>::const_iterator it = stands.constFind(this);
    while (it != stands.constEnd() && it.key()==this) {
        FMStand *stand = it.value();
        stand->reload();
        stand->calculateMAI();
        // calculate sustainable total harvest (following Breymann)
        double area = stand->area();
        mai += stand->meanAnnualIncrementTotal() * area; // m3/yr
        age += stand->absoluteAge() * area;
        volume += stand->volume() * area;
        // HDZ: "haubarer" average increment: timber that is ready for final harvest
        if (stand->readyForFinalHarvest())
            hdz += stand->volume() / stand->absoluteAge() * area; //(0.1* stand->U()) * area; // note: changed!!!! was: volume/age * area
        total_area += area;
        ++it;
    }
    // reset
    ForestManagementEngine::instance()->scriptBridge()->treesObj()->setStand(0);
    mTotalArea = total_area;
    if (total_area==0.)
        return;

    mai /= total_area; // m3/ha*yr area weighted average of annual increment
    age /= total_area; // area weighted mean age
    hdz /= total_area; // =sum(Vol/age * share)

    mMAI = mai;
    //mMAI = mMAI * 1.15; // 15% increase, hack WR
    mHDZ = hdz;
    mMeanAge = age;
    mTotalVolume = volume;

    double rotation_length = U();
    double h_tot = mai * 2.*age / rotation_length;  //
    double h_reg = hdz * 2.*age / rotation_length;
    h_reg = h_tot * 0.85; // hack!
    h_reg *= agent()->schedulerOptions().harvestIntensity;
    h_tot *= agent()->schedulerOptions().harvestIntensity;
    double h_thi = qMax(h_tot - h_reg, 0.);

    qCDebug(abe) << "plan-update for unit" << id() << ": h-tot:" << h_tot << "h_reg:" << h_reg << "h_thi:" << h_thi << "of total volume:" << volume;
    double sf = mAgent->useSustainableHarvest();
    // we do not calculate sustainable harvest levels.
    // do a pure bottom up calculation
    double bottom_up_harvest = (plan_final / period_length) / total_area; // m3/ha*yr

    // the sustainable harvest yield is the current yield and some carry over from the last period
    double sustainable_harvest = h_reg;
//    if (mAnnualHarvestTarget>0.) {
//        double delta = realized/(total_area*period_length) - mAnnualHarvestTarget;
//        // if delta > 0: timber removal was too high -> plan less for the current period, and vice versa.
//        sustainable_harvest -= delta;
//    }
    mAnnualHarvestTarget = sustainable_harvest * sf + bottom_up_harvest * (1.-sf);
    mAnnualHarvestTarget = qMax(mAnnualHarvestTarget, 0.);

    mAnnualThinningTarget = (plan_thinning / period_length) / total_area; // m3/ha*yr


    if (scheduler())
        scheduler()->setHarvestTarget(mAnnualHarvestTarget, mAnnualThinningTarget);
}

QMutex _protect_agent_exec;
void FMUnit::runAgent()
{
    QMutexLocker m(&_protect_agent_exec); // avoid parallel execution of agent-code....

    // we need to set an execution context
    FMStand *stand = *ForestManagementEngine::instance()->stands().find(this);
    if (!stand)
        throw IException("Invalid stand in FMUnit::runAgent");

    FomeScript::setExecutionContext(stand, true); // true: add also agent as 'agent'

    QJSValue val;
    QJSValue agent_type = agent()->type()->jsObject();
    if (agent_type.property("run").isCallable()) {
        val = agent_type.property("run").callWithInstance(agent_type);
        qCDebug(abe) << "running agent-function 'run' for unit" <<  id() << ":" << val.toString();
    } else {
       qCDebug(abe) << "function 'run' is not a valid function of agent-type" << agent()->type()->name();
    }

}

void FMUnit::updatePlanOfCurrentYear()
{
    if (!scheduler())
        return;

    if (mTotalArea==0.)
        throw IException("FMUnit:updatePlan: unit area = 0???");

    // compare the harvests of the last year to the plan:
    double harvests = mRealizedHarvest - mRealizedHarvestLastYear;
    mRealizedHarvestLastYear = mRealizedHarvest;
    mAnnualHarvest = harvests;

    // difference in m3/ha
    double delta = harvests/mTotalArea - mAnnualHarvestTarget;
    mTotalPlanDeviation += delta;

    // apply decay function for deviation
    mTotalPlanDeviation *= mAgent->schedulerOptions().deviationDecayRate;

    // relative deviation: >0: too many harvests
    double rel_deviation = mAnnualHarvestTarget? mTotalPlanDeviation / mAnnualHarvestTarget : 0;

    // the current deviation is reduced to 50% in rebounce_yrs years.
    double rebounce_yrs = mAgent->schedulerOptions().scheduleRebounceDuration;

    double new_harvest = mAnnualHarvestTarget * (1. - rel_deviation/rebounce_yrs);

    // limit to minimum/maximum parameter
    new_harvest = qMax(new_harvest, mAgent->schedulerOptions().minScheduleHarvest);
    new_harvest = qMin(new_harvest, mAgent->schedulerOptions().maxScheduleHarvest);
    scheduler()->setHarvestTarget(new_harvest, mAnnualThinningTarget);

}

} // namesapce