Marie Frapin
email: marie.frapin@mines-paristech.fr

François Chaplais

Patrick Schalbart
email: patrick.schalbart@mines-paristech.fr

Bruno Peuportier
email: bruno.peuportier@mines-paristech.fr

Optimal control of heating in a two-zone building using price decomposition-coordination method

Keywords: Continuous time, Heating management, Indirect method, Optimal control, Two-zone building

. The originality of this paper is to extend this continuous resolution method to a two-zone scale. This is a first step toward multi-zone buildings control. At multi-zone scale, three resolution approaches are possible. The first is the centralised method which consists in solving a unique global OCP. However centralised approach is not suitable for large scale systems due to high computational time or a lack of modularity. At the opposite, the decentralised approach which consists in solving local OCPs at each zone scale is the most used approach nowadays. However, this approach neglects couplings between zones as thermal couplings (type 1) or energy sharing (type 2). The third approach is the decomposition-coordination one which consists in solving local OCPs at each zone scale (decomposition stage) and reintegrate neglected couplings by an iterative procedure (coordination stage). This paper deals with couplings of type 1. The decompositioncoordination method which appeared to be more appropriated was price decomposition-coordination method. It was tested and compared to the optimal reference obtained with the centralised approach in a two-zone building with residential and tertiary zones. The centralised and the decomposition-coordination methods results are close. Both methods leaded to more than 20 % of savings on the global energy bill. The computational time was heavier with the decomposition-coordination method than with the centralised methods. However, in future work concerning multi-zone building, decomposition-coordination method is expected to be more beneficial.

Introduction

In France, buildings represent around 45 % of the energy consumption and 25 % of the greenhouse gases emissions [START_REF]Expérimenter la construction du bâtiment performant de demain[END_REF]. For residential buildings two thirds of the energy consumption is dedicated to heating and for tertiary buildings, heating represents half of the energy consumption [START_REF]Le saviez-vous ? Ademe -Available at <http[END_REF]. One third of heating equipments are electric in both sectors (residential and tertiary) [START_REF]Les chiffres clés du bâtiment[END_REF]. It explains winter electrical peak load which are usually produced by less efficient plants with higher CO 2 emissions (e.g. fuel oil power plant). To avoid electric peak loads, demand side management consists in deferring energy consumption from peak hours to off-peak hours. On/off or PID controllers which are mostly used today for heating management are not efficient to properly decreased peak hours consumption. Indeed, these controllers cannot anticipate time variation of electricity (used to encourage off-peak hours consumption) or internal and external perturbations.

The literature overcomes gaps in classical controllers using advanced regulation methods based on observations as neural network or fuzzy logic [START_REF] Ferreira | Neural networks based predictive control for thermal comfort and energy savings in public buildings[END_REF][START_REF] Marvuglia | Coupling a neural network temperature predictor and a fuzzy logic controller to perform thermal comfort regulation in an office building[END_REF]. Another way is to use Model Predictive Control (MPC) which integrates thermal model of buildings. In recent years, MPC became more famous due to more effective solvers and more accurate physical models especially for heating and cooling systems [START_REF] Ma | Model Predictive Control for the Operation of Building Cooling Systems[END_REF][START_REF] Henze | Experimental Analysis of Model-Based Predictive Optimal Control for Active and Passive Building Thermal Storage Inventory[END_REF]. MPC principle is to solve an optimal control problem (OCP) over a given horizon. Only the first sequence of the optimal solution is apply to the real system and a new optimal control problem is solved with current measurements of the system state and updated data (internal and external loads for instance). For large scale systems, [START_REF] Moroşan | A distributed MPC strategy based on Benders' decomposition applied to multi-source multi-zone temperature regulation[END_REF] shows that computing time increases exponentially with the number of building's zone and promoted distributed MPC. It consists in solving local OCPs (one for each subsystem) and couplings between sub-systems are ensures by a coordination stage. In building energy management, couplings between zones can be thermal couplings (type 1) and energy sharing (type 2). Type 1 were taken into account in [START_REF] Moroşan | Building temperature regulation using a distributed model predictive control[END_REF] and type 2 in [START_REF] Moroşan | A distributed MPC strategy based on Benders' decomposition applied to multi-source multi-zone temperature regulation[END_REF][START_REF] Moroşan | Building temperature regulation using a distributed model predictive control[END_REF][START_REF] Lamoudi | Distributed constrained Model Predictive Control based on bundle method for building energy management[END_REF]. All of these studies were implemented in discrete time. To the authors' knowledge, studies on large scale MPC for buildings energy management using continuous time OCP resolution method are almost inexistent. Continuous time OCP resolution method guarantees the respect of state constraints (comfort levels) at any time which is not always possible with discrete time resolution method especially if time step is large to limit computational times. Continuous time OCP resolution has proved its worth for optimal control of heating in a monozone building using on indirect resolution method based on Pontryagin's minimum principle [START_REF] Robillart | Etude de stratégies de gestion en temps réel pour des bâtiments énergétiquement performants[END_REF]. The aim of this paper is to apply price decomposition-coordination method to optimal control of heating in a two-zones building (a residential zone and an office zone) using a continuous time method. Results are compared to the centralised optimum. The paper is organised as follows. Section 2 introduces the OCP statement. Sections 3 and 4 respectively present the centralised and decomposition-coordination methods to solve the OCP. In section 5, both methods are applied and discuss on a study case. Section 6 concludes the study and outlines future works.

Optimal control problem statement

The optimal control principle is to find a control law over a given time horizon such that an optimal objective is achieved. Following the presentation of the objective function for a two-zones building, all the considered constraints of the OCP are depicted in this section.

Objective function

The overall control objective is to minimise the heating cost of a two-zones building: [START_REF] Robillart | Etude de stratégies de gestion en temps réel pour des bâtiments énergétiquement performants[END_REF] where u is the controls vector containing the electric heating power of each zone, represents the time dependent electricity cost, and are respectively the initial and the final times.

Dynamic system modelling

For modelling sake, a building with an arbitrary number of zones is broken down into a set of meshes (finite volumes). A thermal balance on each mesh element is applied to establish the overall dynamic system defined as follows: [START_REF]Expérimenter la construction du bâtiment performant de demain[END_REF] where x is the state vector containing the temperature of each mesh, s is the vector of internal and external loads and y is the outputs vector containing the indoor temperature of each zone. Matrices A, B, C and D are adequately sized.

The equations system (2) represents a dynamic constraint of the OCP which ensures the respect of the building thermal behaviour.

State constraints

To provide comfort for occupants, state constraints are considered to ensure that indoor temperature of each zone i ( ) is within a fixed interval:

(3) Equation ( 3) can be rewritten depending on the state x:

(4)
The state constraints define the OCP as a dynamic problem. Indeed, without these, no need to heat the building.

Control constraints

Due to production systems limitations, constraints on control were considered to ensure that the heating power of each zone i ( ) is also within a fixed interval characterised by:

(5)

Centralised method

After defining the OCP, the focus was on solving it in continuous time. Before setting up the decomposition-coordination resolution method, the centralised method was implemented to find the reference solution. The centralised method is depicted in this section, it helps introducing the following section about decomposition-coordination method. The centralised method is based on the algorithm developed by [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of nonlinear systems[END_REF] for monozone buildings' OCP. This algorithm handles state and control constraints using internal penalty method described in the first part of this section. In this work, state constraints were effectively implemented by internal penalty method. However, control constraints were directly integrated using Pontryagin's minimum principle (second part of this section). The centralized algorithm is described in the last part of this section.

Internal penalty method for state constraints implementation

The principle of internal penalty methods consists in penalising the cost function by a barrier function which strongly increased the cost when the solution approaches the state constraint. With an example of state constraints defined as: [START_REF] Marvuglia | Coupling a neural network temperature predictor and a fuzzy logic controller to perform thermal comfort regulation in an office building[END_REF] the penalised cost function takes the following form: [START_REF] Ma | Model Predictive Control for the Operation of Building Cooling Systems[END_REF] with typical defined as: [START_REF] Henze | Experimental Analysis of Model-Based Predictive Optimal Control for Active and Passive Building Thermal Storage Inventory[END_REF] The penalised problem is now state constraints free. Its solution converges to the optimal one with state constraint as tends to 0. Considering that the centralised problem has four state constraints of the form defined by inequality (6) (two for each zone), the penalised cost function of the centralised two-zones problem is:

(9)

Optimal control research with control constraints (minimisation of the Hamiltonian)

According to Pontryagin's minimum principle (Appendix A), the optimal triplet must satisfy the following assertion: [START_REF] Moroşan | Building temperature regulation using a distributed model predictive control[END_REF] where p is the costate satisfying: [START_REF] Lamoudi | Distributed constrained Model Predictive Control based on bundle method for building energy management[END_REF] and the Hamiltonian defined by: [START_REF] Mesarovic | Two coordination principles and their application in large scale systems control -ScienceDirect[END_REF] To simplify Hamiltonian's writing, time does not appear in equation [START_REF] Mesarovic | Two coordination principles and their application in large scale systems control -ScienceDirect[END_REF]. To easily solve [START_REF] Moroşan | Building temperature regulation using a distributed model predictive control[END_REF], a small quadratic term in u is added to the objective function which leads to the new Hamiltonian: [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of nonlinear systems[END_REF] where is a positive real sized to ensure that the quadratic term does not influence the heating cost: [START_REF] Carpentier | Décomposition-coordination en optimisation déterministe et stochastique[END_REF] The Hamiltonian defined in equation ( 13) is a polynomial of degree two in u. The solution of equation ( 10) is either the parabola vertex or one of the limits of interval .

Centralised algorithm

The centralised resolution algorithm includes the three following steps.

 Step 1: initialisation

The state function is initialised in the internal set of state constraints (4) for all and the costate function is initialised to 0 for all . The decreasing sequence is defined and initialized at .

 Step 2: two boundaries problem resolution with inner minimisation of the Hamiltonian

The optimal control is requested minimising the Hamiltonian with respect to u.

The two boundaries problem consists in solving the following equations system:

(15)

 Step 3: increment

If the optimum is reached and the algorithm stop, else is incremented and the algorithm returns in step 2 with the new values of and as initialisation.

Price decomposition-coordination method

The second approach of global OCP solving discussed in this paper is a decomposition-coordination method. Decomposition-coordination method consist in solving local OCPs (at zone scale) in a decomposition stage and ensuring couplings between zone via iterations with a coordinator in a coordination stage. Thermal couplings between zone lead to couplings between dynamical systems of each zone. [START_REF] Mesarovic | Two coordination principles and their application in large scale systems control -ScienceDirect[END_REF] firstly introduced the interaction balance principle to deal with coupled dynamic systems. [START_REF] Carpentier | Décomposition-coordination en optimisation déterministe et stochastique[END_REF] proposed three resolution methods in is second model dedicated to coupled dynamic systems (the first model is rather dedicated to energy sharing case). Among this resolution methods, there is resource allocation method (or primal method) and the price method (or dual method).

While the resource-allocation approach ensures the strict respect of couplings at each iteration, we focused on price method where in comparison to primal method, no gradient need to be provided in local optimization (decomposition stage). Indeed, computations of gradient is CPU intensive. The third method called prediction method can mixes both primal and dual method. It was not selected for the same reason as for the primal method.

Optimal control sub-problem statement

System (2) can be decomposed into two sub-systems associated to each zone i and defined by the following dynamic system:

(16)

The purpose of decomposition is to solve independent optimal control sub-problem in order to parallelise the resolution. As mentioned above, sub-systems are dynamically coupled. To make subsystems independent, for each zone, the term is replaced by an interaction variable :

(17)
The interaction variable becomes a new control of zone i .Then the coordination purpose is to ensure coupling constraints defined as follows:

(18)

To do that, prices and are added to the global objective function ( 1) in order to force the solution to respect constraint defined by equation ( 18). The modified global Hamiltonian is defined as follows: For simplicity sake, state constraints are not specified (hence the superscript 0). They are identical to the one in equation [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of nonlinear systems[END_REF]. A saddle point is supposed so the optimum must verified the following assertion:

(20)

Particularly, the minimisation problem can be rewritten as: Price decomposition-coordination method is based on Uzawa algorithm which is in fact a gradient algorithm on multipliers and at the coordination stage [START_REF] Carpentier | Décomposition-coordination en optimisation déterministe et stochastique[END_REF].

Augmented Lagrangian

Uzawa algorithm is facilitated in linear case using augmented Lagrangian. It consists in adding the following quadratic term to the objective function : (23) At iteration k of the decomposition-coordination algorithm, the augmented Lagrangian additional term ( 23) is divided into two components assigned to each local Hamiltonian: (24) with k the decomposition-coordination index.

Decomposition stage

Decomposition stage goal is to minimise local Hamiltonians in a parallel way. Each local OCP is controlled in and .

is constrained by the fact that is constraineed. The variation range of is noted and is determined by equations ( 4) and (18). According to Pontryagin's minimum principle, the optimal set must satisfy the following assertion:

Local Hamiltonians are polynomials of degree two in and in . The solution of equation ( 25) is either the parabola vertex in or one of the limits of interval . The solution of equation ( 25) is either the parabola vertex in or one of the limits of interval .

Coordination stage

Coordination stage goal is to update as follows:

(

Convergence is ensured when the following assertion is checked:

(27)

Decomposition-coordination algorithm

 Level 1: global initialisation Prices and are initialised for all . We fixed the index k at 0.

 Level 2: decomposition

For each zone, the three following steps are done. The three steps are similar to those defined in centralised algorithm except that controls are now and .

 Step 1: local initialisation

Functions is initialised to respect all the state constraints defined by inequality (4) for all and function is initialised to 0 for all . The decreasing sequence suits is defined and initialised at .

 Step 2: two boundaries problem resolution with inner minimisation of the Hamiltonian

The optimal controls and are requested to minimise the Hamiltonian with respect to and . The two boundaries problem consists in solving the following equations system:

Application of the optimal control for both methods

Case study

The case study was a two-zones building, one for residential use (zone 1) and one for offices (zone 2). Both zones were considered exposed to the same external loads (Fig. 2a) but was characterised by different occupation times (from 5 p.m. to 9 a.m. for zone 1 and from 9 a.m. to 7 p.m. for zone 2). This led to the internal heat gains profiles presented at Fig. 2.b.

Fig. 2. a) External loads, b) internal loads

States constraints ensured minimal comfort temperature setpoints of 21 °C for zone 1 and 23 °C for zone 2. Outside of occupation times, reduced temperature setpoints equal to 18 °C for zone 1 and 20 °C for zone 2 were taken into account. Maximal temperature setpoints equal to 26 °C for zone 1 and 28°C for zone 2 were also considered. Regarding control constraints, for each zone, was constrained to be between 0 and 20 kW. Three electricity cost were considered depending on the hour. Between 0 a.m. and 9 a.m. the electricity price was 0.1 €/kWh, this period is called off-peak hours. This price was doubled between 9 a.m and 5 a.m and between 10 p.m and 0 a.m., this period is called peak hours. The initial price was tripled during high peak hours between 5 p.m and 10 p.m.

Implementation

Computations were performed on a 8-Core processor at 2.80 GHz, with 32 GO RAM and running Windows 7 Professional. Both centralised and decomposition-coordination methods were implemented in Matlab using bvp5c routine to solve two-boundaries problems. In decomposition stage of decomposition-coordination method, the computations were parallelised with the parallel computing toolbox.

Algorithm parameters

Time horizon was 3 days, however, only the first two days are represented in the results. Indeed, a boundary-effect was observed on the last day of optimisation, induced by the impossibility on the last day to anticipate the needed storage effort for the following days. For both centralised and decomposition-coordination methods, the sequence was defined to make varying between 10 -5 and 10 -7 . For the decomposition-coordination method, the tolerance tol mentioned in equation ( 27) was equal to 0.15 °C. Experiments showed that a good estimation of the coordination step and the augmented Lagrangian parameter r was 2.10 -8 . Regarding prices and , they were initialised at 2.10 -5 for all .

Centralised method results

Fig. 3 represents the optimal heating power profiles (a) and the resulting temperature profiles (b) of both zones 1 and 2 achieved with the centralised method. Each zone is pre-heated during off-peak hours to shave off heat consumption during high peak hours, also ensuring that temperature setpoints are respected.

In comparison with the case without demand side management, the total cost was reduced by 22.9 %. For zone 1, 99 % of energy consumption during high peak hours and 100 % of energy consumption during peak hours were reduced. For zone 2, 17 % of energy consumption during high peak hours and 57 % of energy consumption during peak hours were reduced. The minimisation of peak loads is more important for zone 1 (residential space) than for zone 2 (office space). Indeed, the periods when zone 1 needs further heat correspond essentially to off-peak hours unlike zone 2, where further heat needs are exclusively during peak and high peak hours.



  Step 3: increment If the optimum is reached and the algorithm stop, else is incremented and the algorithm returns in step 2 with the new values of and as initialisation  Level 3: coordination Prices and are updated with equation (26). Level 4: global incrementWhile condition (27) is not achieved, increment returned to step 2. Fig.1illustrates the decomposition-coordination algorithm for better understanding.

Fig. 1 .

 1 Fig. 1. Decomposition-coordination algorithm

Fig. 3 .

 3 Fig. 3. Centralised method results : (a) optimal heating power profiles, (b) resulting temperature profiles.

Price decomposition-coordination method results

The three decomposition-coordination methods detailed in [START_REF] Carpentier | Décomposition-coordination en optimisation déterministe et stochastique[END_REF] were implemented in this study. As mentioned above, computations for resource allocation and prediction were too heavy due to the need to provide gradients in local OCPs resolution. The price decomposition-coordination method was preferred. Its results are presented in this section. Fig. 4 represents the optimal heating power profiles (a) and the resulting temperature profiles (b). While local OCPs resolution are done with local dynamic systems (17), resulting temperature profiles were obtained by feeding back optimal heating power profiles into the complete model [START_REF]Expérimenter la construction du bâtiment performant de demain[END_REF]. This procedure allows to confirm that comfort levels are respected. Indeed, Fig. 4 (b) shows that comfort levels are respected even if temperature profiles get closer to the temperature setpoints compared with the centralised method results (Fig. 3 (b)). In price decomposition-coordination algorithm, the cost was reduced by 22,6 %. For zone 1, 97 % of energy consumption during high peak hours and 100 % of energy consumption during peak hours were reduced. For zone 2, 15 % of energy consumption during high peak hours and 55 % of energy consumption during peak hours were reduced. As for centralised method results, the minimisation of peak loads is more important for zone 1 than for zone 2. Regarding cost and energy reduction, centralised method results are slightly better than decomposition-coordination method results.

Fig. 4. Decomposition-coordination method results : (a) optimal heating power profiles, (b) resulting temperature profiles

Conclusion and future work

To conclude, two continuous time resolution methods were performed to find the optimal control of heating in a two-zone building mixing residential and office spaces. The first method applied was the centralised method. This method which consists in solving a unique global OCP converges to the global optimum. To anticipated the increasing complexity of the centralised method at multizone scale, a second method called price decomposition-coordination method was applied. Results of the decomposition-coordination method were compared to the global optimum. Regarding cost and energy reduction, results are close even if they are slightly better with the centralised method. With both methods, high peak and peak hours consumptions declined by almost 100 % for zone 1. For zone 2, high peak hours consumption were reduced almost in half nevertheless, peak hours consumption reduction did not exceed 20%. The decomposition-coordination method presented higher computational time than centralised method. In the experiments, price and were initialised to a constant value. Future work consists in implementing an MPC loop. In this case, prices profiles and can be reused from one OCP resolution to the next one. This is expected to facilitate the decompositioncoordination algorithm convergence and therefore to reduce computational time. Finally the generalisation of both methods to multi-zone buildings is envisaged . Indeed, while computations of centralised method increased exponentially with the number of zones, decomposition-coordination method computational time is expected to increase less dramatically thanks to parallelisation.

Appendix A

This Appendix details the Pontryagin's minimum principle. This principle is a necessary optimality condition asserting that if the control associated to the state is the solution of the OCP then it exists an application called adjoint state which fulfil following assertions:

This necessary optimality condition leads to solve the following system of equations:


Emmanuel Paradis
email: paradis@isem.univ-montp2.fr

Julien Claude

Korbinian Strimmer

Keywords: R, computer language, molecular evolution, phylogenetics, statistics

APE (Analysis of Phylogenetics and Evolution) is a package written in the R language for use in molecular evolution and phylogenetics. APE provides both utility functions for reading and writing data and manipulating phylogenetic trees, as well as several advanced methods for phylogenetic and evolutionary analysis (e.g. comparative and population genetic methods). APE takes advantage of the many R functions for statistics and graphics and also provides a flexible framework for developing and implementing further statistical methods for the analysis of evolutionary processes.

Phylogenetic analysis, in its broad sense, covers a very wide range of methods from computing evolutionary distances, reconstructing gene trees, estimating divergence dates, to the analysis of comparative data, estimation of evolutionary rates, and analysis of diversification. All of these diverse tasks have one particular aspect in common: they rely heavily on computational statistics.

The R system, a free platform-independent open-source analysis environment, has recently emerged as the de-facto standard for statistical computing and graphics [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF]. One advantage of R is that it can be easily tailored to a particular application area by writing specialized packages. In particular, the usefulness of R in bioinformatics has already been impressively demonstrated in the analysis of gene expression data (http://www.bioconductor.org).

APE (Analysis of Phylogenetics and Evolution

) is the first joint effort to utilize the power of R also in the analysis of phylogenetic and evolutionary data. APE focuses on statistical analyses using phylogenetic and genealogical trees as input. In version 1.1, APE provides functions for reading, writing, plotting, and manipulating phylogenetic trees, analyses of comparative data in a phylogenetic framework, analysis of diversification, computing distances from allelic and nucleotide data, reading nucleotide sequences, and several other tools such as Mantel's test, computation of minimum spanning tree, or estimation of population genetics parameters. Table 1 gives an overview of the functions currently available in APE.

Note that some of the methods (e.g., comparative method, skyline plot, etc.) have previously been available only in specialized softwares. External tree reconstruction programs (such as PHYLIP) can be called from R through standard shell commands.

One strength of R is that it is straightforward to obtain publication-quality graphical output, particularly with its PostScript device. For instance, the plotting function of phylogenies in APE handles colors, line thickness, font, spacing of labels, which can be defined separately for each branch, so that three different variables can be represented on a single phylogeny plot. APE also produces complex population genetics plots, such as the generalized skyline plot [START_REF] Strimmer | Exploring the demographic history of a sample of DNA sequences using the generalized skyline plot[END_REF], with a single command. APE, like any R package, is command-line driven. The functions are called by the user, possibly with arguments and options. Any session using APE in R starts with the command library(ape) which makes the functions of APE available in the R environment. The list of these functions can be displayed with the command library(help = ape) which displays their names with a brief description. An evolutionary tree saved on the disk in the text file tree1.txt in the standard Newick parenthetic format can then be read by tree1 <-read.tree ("tree1.txt") This stores the phylogenetic tree is in an object named tree1 of class "phylo".

The information stored in this object (e.g., branch lengths) can be inspected by typing tree1 and graphical output in form of a cladogram can be obtained by executing plot(tree1) which actually calls the function plot.phylo of APE to draw the phylogenetic tree tree1 (due the object-oriented nature of R the command plot(x) may give a completely different result depending on the class of x). The tree is plotted, by default, on a graphical window, but can be exported in various file formats depending on the operating system.

In addition to this trivial example, the representation of a phylogenetic tree in an object-oriented structure results in straightforward manipulation of the phylogenetic data for various computations used in evolutionary analyses. Currently implemented in APE are approaches such as phylogenetically independent contrasts [START_REF] Felsenstein | Phylogenies and the comparative methods[END_REF][START_REF] Harvey | The Comparative Method in Evolutionary Biology[END_REF], fitting birth-death models [START_REF] Nee | The reconstructed evolutionary process[END_REF][START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF], population-genetic analysis [START_REF] Nee | Inferring population history from molecular phylogenies[END_REF][START_REF] Strimmer | Exploring the demographic history of a sample of DNA sequences using the generalized skyline plot[END_REF], non-parametric smoothing of evolutionary rates [START_REF] Sanderson | A nonparametric approach to estimating divergence times in the absence of rate constancy[END_REF], and estimation of groups of genes in phylogenetic trees using Klastorin's method [START_REF] Misawa | A simple method for classifying genes and a bootstrap test for classifications[END_REF]. Furthermore, distance-based clustering methods as implemented in the R function hclust can be used by APE using functions converting to and from objects of class "phylo" and "hclust".

All R functions available in APE (see Table 1) are documented in the R hypertext format and information regarding their use can be retrieved by applying the help command, e.g., help(read.tree)

The classes and methods in APE (like phylo) can also easily be further extended to include other functionalities, for instance to annotate phylogenetic trees. Thus, APE is not only a data analysis package, it is also an environment for developing and implementing new methods. Furthermore, programs written in C, C++, or Fortran77 can be linked and called from R. This is particularly useful for computer intensive calculations.

Table 1 . Special functions available in APE 1.1

 1 Miscellaneousall.equal.phylo, balance, base.freq, dist.dna, dist.gene, dist.phylo, GC.content, klastorin, mantel.test, mst, summary.phyloDetailed information about each function can be accessed with the online help (e.g.

	Application	Available commands
	Input/Output	read.dna, write.dna, read.nexus,
		write.nexus, read.tree, write.tree,
		read.GenBank
	Graphics	add.scale.bar, plot.mst, plot.phylo,
		plot.skyline, lines.skyline, ltt.plot
	Tree manipulation	bind.tree, drop.tip,
		is.binary.tree, is.ultrametric
	Comparative method compar.gee, compar.lynch, pic, vcv.phylo
	Diversification	birthdeath, cherry, diversi.gof,
		diversi.time, gamma.stat
	Population genetics	branching.times, coalescent.intervals,
		collapsed.intervals, find.skyline.epsilon,
		heterozygosity, skylineplot,
		skyline, theta.h, theta.k, theta.s
	Molecular dating	chronogram, ratogram, NPRS.criterion
	Data sets	bird.families, bird.orders, hivtree,
		landplants, opsin, woodmouse, xenarthra

help(mantel.test)).

Acknowledgments

We thank two anonymous referees for their constructive comments on a previous version of this paper. This research was financially supported by the Programme inter-EPST "Bioinformatique" (E.P. and J.C.) and by an Emmy-Noether research grant from the DFG (K.S.).

The program is free and available from the official R package archive at http://cran.r-project.org/src/contrib/PACKAGES.html#ape.